Aerodynamic-driven topology optimization of compliant airfoils

A strategy for density-based topology optimization of fluid-structure interaction problems is proposed that deals with some shortcomings associated to non stiffness-based design. The goal is to improve the passive aerodynamic shape adaptation of highly compliant airfoils at multiple operating points. A two-step solution process is proposed that decouples global aeroelastic performance goals from the search of a solid-void topology on the structure. In the first step, a reference fully coupled fluid-structure problem is solved without explicitly penalizing non-discreteness in the resulting topology. A regularization step is then performed that solves an inverse design problem, akin to those in compliant mechanism design, which produces a discrete-topology structure with the same response to the fluid loads. Simulations are carried out with the multi-physics suite SU2, which includes Reynolds-averaged Navier-Stokes modeling of the fluid and hyper-elastic material behavior of the geometrically nonlinear structure. Gradient-based optimization is used with the exterior penalty method and a large-scale quasi-Newton unconstrained optimizer. Coupled aerostructural sensitivities are obtained via an algorithmic differentiation based coupled discrete adjoint solver. Numerical examples on a compliant aerofoil with performance objectives at two Mach numbers are presented.


Introduction
Topology optimization represents a radical departure from conventional sizing methods as it allows an optimum material distribution to be identified. It has been applied to aircraft structures for over two decades (Balabanov and Haftka 1996). In most applications, the technique is applied locally, e.g., to single ribs (Krog et al. 2004), so that the resulting structure can still be manufactured by traditional methods. Another important practical challenge of topology optimization, especially in a fluid-structure interaction (FSI) context, is the computational cost. While sometimes the analysis can be simplified by assuming "frozen" fluid loads (see Zhu et al. 2016), this assumption can lead approach, e.g., using the SIMP (solid isotropic material with penalization, Bendsøe 1989) or RAMP (rational approximation of material properties, Stolpe and Svanberg 2001) interpolation schemes, as an optimum solid-void topology is contingent on how important the local stiffness to weight ratio of the ersatz material is to the objectives. The material interpolation is such that intermediate density areas have lower stiffness to weight ratios than those of solid or void areas (the latter have some residual stiffness); therefore, such areas are undesirable in an optimally stiff structure. In general, aerodynamic objectives do not have this property and so a careful formulation that recovers it may be required, and explicit penalization of nondiscreteness has also been proposed (Stanford and Ifju 2009). Alternative topology optimization methods such as the level-set approach (e.g., Dunning et al. 2015) inherently produce a solid-void topology. However, the numerical challenges of the density-based approach are better understood (e.g., control of feature sizes) and their natural ability to introduce new holes, and be driven by a general optimizer (which easily allows other variables to be included in the optimization) is not shared by all level-set approaches.
In this work, we present a method to apply density-based topology optimization to a highly compliant airfoil, with the objective of improving its aerodynamic characteristics at multiple operating points. This differs from most previous topology optimization research in that the objectives are mostly aerodynamic, which for the reasons above, requires high-fidelity modeling of the fluid to be used. Furthermore, hyper-elasticity is considered on the structural side as strains are large and so accounting for nonlinear behavior (including buckling) may be necessary to accurately model the structural response. Due to the nature of the objectives and the large number of design variables, gradient-based optimization is used. The computational platform used is the multi-physics suite SU2 ) which has been verified for a range of aeronautical applications, e.g., in Palacios et al. (2014). The FSI primal and adjoint solution methodology is presented in Section 2 and the optimization strategy in Section 3. The numerical results obtained with the proposed strategy are presented in Section 4, and finally in Section 5, they are compared with results obtained using common ways to encourage solid-void solutions in densitybased methods, and we discuss why those failed due to the characteristics of the example problem.

Primal and adjoint coupled solution methods
A three-field partitioned formulation is adopted for the FSI problems (primal and adjoint). Coupled sensitivities are obtained with the discrete adjoint method. This is outlined next together with some details of the computational implementation in SU2.

Fluid dynamics
In the fluid domain, the flow is governed by the continuity, Navier-Stokes, and energy conservation equations, which in the ALE formulation may be written as where w = (ρ, ρv, ρE) is the vector of conservative variables. The convective and diffusive fluxes, and the volumetric sources are defined, respectively, as where ρ is the density, v andż, the flow and grid velocities in a Cartesian coordinate system, respectively (z are the grid displacements), p the pressure, E the total energy per unit mass, τ the stress tensor, κ the thermal conductivity, and T the temperature. All material properties refer to the fluid. For the viscous stress tensor τ , Newtonian fluid is assumed and bulk viscosity effects are ignored; furthermore, under the Boussinesq hypothesis, turbulence is modeled as increased viscosity. Menter's Shear Stress Transport turbulence model (Menter 1993) is used. Pressure and temperature are related with the conservative variables via the ideal gas equation of state.
Equation (1) is integrated in space using the finite volume method (FVM) on a dual grid with control volumes constructed using a median-dual vertex-based scheme, which results in the following semi-discrete integral equation for each volume i where the residual R i results from summing the discretized fluxes for all faces of the control volume and integrating the volumetric sources; we use the JST convective scheme for its robustness knowing that the values of drag coefficient will be overestimated due to artificial dissipation. To obtain a steady-state solution, (5) is marched implicitly in pseudo time (τ ), that is, the new solution w * is obtained by solving where the continuous temporal derivative has been replaced by a backward-Euler approximation and the tilde indicates that the linearization of the residual is approximate. For example, more weight is given to the Jacobian of the dissipation term to mitigate the non-diagonal dominance that results from central schemes.
We represent the fluid solution process, involving at each iteration the computation of the residual and the solution of the linear system in (6), by the fixed-point iteration The turbulence model equations are solved in the same manner but they are lagged; nevertheless, for the purposes of the fixed-point representation, the turbulence variables can be considered part of w.

Solid mechanics
The deformed state of a solid domain ( ) is governed by the point-wise equilibrium of linear momentum and of tractions on its surface ( ), that is, where ρ is the density,ü the acceleration with respect to an inertial frame, f the inertial body forces, σ the Cauchy stress tensor, n the outward surface unit normal, and λ the external tractions. To solve (8) for the structural displacements via the finite element method (FEM), its weak form is established applying the principle of virtual work (see Bonet and Wood 2008), where δE is the variation of the Green-Lagrange strain tensor with respect to δu, S is the second Piola-Kirchhoff stress tensor, r refers to the reference (undeformed) configuration of the structure, and c and c to the current configuration of the structure and its surface, respectively. For hyper-elastic materials, the relation between stress and strain is given by the strain energy density function , as S = ∂ ∂E . In particular, for a neo-Hookean material with Lamé constants μ and λ, it is given as where J is the determinant of the deformation gradient.
Equation (10) is linearized around the current state of deformation before being discretized, linear isoparametric elements are used in SU2 (Sanchez et al. 2016). The solution is then found iteratively via the Newton-Raphson method which again can be formulated as a fixed-point iteration u = S(u, λ(w, z)) (12)

Mesh deformation
The structural displacements at the fluid-structure interface (u ) are transferred to the interior nodes of the fluid domain using a linear elasticity analogy (Sanchez et al. 2016). This can be stated as the explicit relation (13), which represents the third field of the FSI problem.

Load and displacement transfer
In general, the fluid and structural meshes do not match at their interface, and an interpolation scheme is required to transfer displacements from the solid to the fluid side, and tractions in the opposite direction. A radial basis function (RBF) approach (see Beckert and Wendland 2001) is adopted in this work to generate the required interpolation matrices (H), that is where subscripts s and f stand for solid and fluid, respectively, and the fluid tractions are given by with n f the inward (with respect to the fluid domain) surface unit normal vector.

Partitioned coupling method
Under a suitable level of convergence of the fixedpoint iterations, the fluid solver can be considered the Dirichlet-to-Neumann operator F that maps displacements to tractions, and the structural solver the Neumann-to-Dirichlet operator S that does the opposite, i.e., λ f = F(u ) and u = S(λ f ) where for simplicity of notation the transfer and mesh deformation operations have been lumped with the solvers. These operators can be applied sequentially to produce a fixed-point iteration for the interface displacements, which naturally results in a block-Gauss-Seidel (BGS) approach to solve the interface problem. However, this often requires significant under-relaxation for strongly coupled problems. Thus, the interface quasi-Newton method IQN-ILS (Degroote et al. 2009) was also implemented in SU2 for this work, to accelerate convergence and improve the robustness of the coupled solver. The interface residual is given by and the new interface position, u * , by Note that u * =û , even for BGS if some relaxation is applied. The product of the (least-squares) approximate Jacobian inverse (r −1 u ) with the residual is obtained based on (up to) N previous values of r andû , that is where each column of W stores the difference between consecutive values ofû , that is and the coefficients c are obtained by solving the linear least-squares problem where analogously to W, V = r N − r N−1 , · · · , r 2 − r 1 , and the yet unknown residual at the next coupling iteration (r N+1 ) is desired to vanish and therefore assumed to be 0. We have observed that for the problems in this work, IQN-ILS is 1.5 times faster than BGS.

Coupled sensitivities
The coupled adjoint solution algorithm follows on the work of Albring et al. (2016) and Sanchez et al. (2018). Let G represent the fixed-point iterator for the coupled problem consisting of the concatenation of (7), (12), and (13), and let x = (w, u, z) represent the state of the coupled problem and α the parameters with respect to which some functional J is to be minimized, i.e., introducing the adjoint variables x = (w, u, z) the Lagrangian is defined as Differentiating with respect to the parameters and grouping some terms Defining the adjoint variables such that the term in parenthesis vanishes results in the adjoint fixed-point iterator which, considering the 3-field nature of the problem and the direct dependencies of each operator, becomes ⎛ where subscripts indicate differentiation with respect to the given variable (e.g., F w ≡ ∂ w F). Analogous to what was done for the primal problem, iterating on each adjoint variable with the remaining fixed allows one to define adjoint solvers for each field, that is Noting that by the construction of M, it is M u = 0 ∀ u / ∈ ; the adjoint interface displacements are identified as which allows defining the adjoint interface residual as The IQN-ILS method can then be applied to determine u and therefore x, which upon substitution into (25) yields the sensitivities. In SU2, the linearization of the fixedpoint operators around the converged state, to compute matrix-free Jacobian-transposed products, is done using the algorithmic differentiation (AD) tool CoDiPack (Sagebaum et al. 2019).

Optimization strategy
As mentioned in the introduction, aerodynamic optimization objectives generally do not benefit from an optimum stiffness to mass ratio, and therefore intermediate densities may not be avoided using the SIMP formulation. This is especially true when seeking passive load alleviation, as deformation is necessary, and in 2D where an increase in lift due to higher mass does not penalize drag as severely as it does in 3D due to higher induced drag. Therefore, an explicit mass objective or constraint may be required. However, we found this to make the optimization process less robust, as a strong-enough incentive to remove material often conduces to critically stable structures (for which it is difficult to fully converge the primal and adjoint problems). Furthermore, we also observed these strategies to increase the sensitivity to initial conditions. For example, starting from an initial configuration that produces more lift than required will cause material to be quickly removed from high strain energy areas. Moreover, mass is not a primary objective for the purposes of this study as we wish to assess the potential of manipulating aerodynamic performance via a material distribution. Therefore, it is not trivial to determine to what value mass should be constrained or how it should be weighted before combining it with the aerodynamic objective function, and numerical examples of this are included in Section 5. To avoid the aforementioned issues, we investigate a two-step process that reduces the interference between the goals of improving aerodynamic performance and producing a realizable topology.

Density-based topology optimization
The well-established density approach with continuous variables consists in specifying a design density at discrete locations of the solid domain, typically the element centroids, and making the local elasticity modulus a function of it. In this work, the modified SIMP formulation (31) is used to relate the elasticity modulus with the design variable as where E min is introduced to avoid a singular stiffness matrix and is typically at least three orders of magnitude smaller than the reference value for solid material (E 0 ), despite this relaxation, a direct sparse solver (Hénon et al. 2002) is required due to the non-Cartesian grids used in the numerical examples. For values of the penalization exponent p greater than 1, intermediate density areas have unfavorable stiffness to weight ratio, and so under the right objectives, constraints, and penalization, will tend to be eliminated. However, for intermediate densities to be realizable (via a two-phase micro-structured material), their properties must respect the Hashin-Shtrikman bounds for two-phase materials (Bendsøe and Sigmund 2004), which in 2D requires p ≥ 3 for the Poisson ratio of 1/3. Two well-known numerical difficulties associated with this approach are its lack of grid convergence (as the discretization is refined, more holes can be introduced) and the checker-boarding that may result from the over estimation of the stiffness of corner contacts. To avoid these issues, a discrete filtering operation (Bruns and Tortorelli 2001) is introduced between the design density variables exposed to the optimizer, and the physical densities considered for each finite element, i.e., where w ij = (R − ||x i − x j ||) j and N is the set of elements within radius R of element i. This conical filter kernel invariably results in intermediate density halo regions between solid and void regions. Filters built by combining the morphological operators (Sigmund 2007) dilatẽ and its counterpart erode (ρ i = 1 − f i (1 − ρ)) were able to produce discrete topologies (at β ≈ 200) for the canonical problems used to verify the implementation. However, numerical investigations for the FSI problem in this work showed they complicate the convergence of the optimization due to the non-linearity they introduce. In either case, derivatives with respect to the design densities are obtained in the adjoint post-processing step (25) using AD (α ≡ ρ); thus, analytical expressions for the Jacobian matrices G α and J α are not required.

Gradient-based optimization method
The numerical optimization examples considered in this work are characterized by a large number of design variables and relatively few constraints (excluding simple bound constraints) that from the physics point of view, do not need to be imposed strictly. Therefore, we solve the optimization problems with the exterior penalty method, i.e., problems of the form where h + = max(0, h). The L-BFGS-B implementation available with SciPy (Jones et al. 2001) is used as the unconstrained (but bounded) optimizer. The penalty parameters (a i and b j ) need to be gradually increased (usually by multiplying the previous value by a fixed factor r) until a predetermined small constraint tolerance is met. This creates the need for outer iterations as updating the parameters within L-BFGS-B (inner) iterations leads to bad approximations of the Hessian matrix. Although these outer iterations force an undesired reversion to steepest descent, they are also needed (and commonly used) to update topology optimization parameters. Our continuation strategy for these parameters, which we have found to be adequate for both structural and FSI problems, is given in algorithm 1.
Note that penalty parameters are increased geometrically, whereas problem parameters are simply a sequence of values through which the algorithm advances if all constraints are currently satisfied. Before the target values of the parameters are reached, loose convergence criteria are used for L-BFGS-B (e.g., 40 inner iterations). The objective function is shifted and scaled by representative minimum value and range, respectively; the constraints are shifted by their bounds and scaled by a reference value (the reciprocal of the bound unless otherwise specified). Doing so allows the same constraint tolerance (≈ 0.01) to be used, the penalty parameters to be initialized equally (a 0 i , b 0 i ∈ [1, 10]), and also updated with the same factor (r ∈ [1.4, 4]).

Strategy for discrete topologies
The proposed two-step strategy consists of first solving the natural aerostructural optimization problem, for example minimizing drag subject to a lift constraint, without additional penalties or manipulation of the objectives. This first step is conducted with fully coupled FSI modeling (so called multidisciplinary-feasible approach) and generally results in a non-discrete topology, but one defined by feasible FSI solutions and realizable since the SIMP exponent is still set according to the two-phase material bounds. The second step then aims to produce a discrete topology that replicates the response of the former, i.e., one that under the fluid loads known from the coupled FSI simulation results in the same deformed surface. This inverse design step is formulated similarly to a compliant mechanism design problem, for example the force inverter (see Bendsøe andSigmund 2004 or Sigmund 2007), but instead of focusing on the response of a single node, an error metric is defined for the entire interface as and used as a constraint in a mass minimization problem.
We solve the problems of both steps with algorithm 1. While it would be possible to perform the inverse design step simulating only the structure, we found that doing so can easily result in an unstable structure that buckles under the small variation of FSI loads caused by the small discrepancy between the target response (obtained from the fully coupled step) and the response of the discrete-topology structure. To mitigate this issue, a one-way FSI simulation is used instead, then based on the variation of the fluid loads due to the discrepancy in target and current surface coordinates, we define a stability metric based on points where the variation of work is positive, that is The one-way simulation is relatively inexpensive as good initial conditions for the flow field are available from the coupled FSI simulation, especially if the initial topology is almost feasible in error metric (35) (the topology obtained in the fully coupled step is not used as the starting point to avoid convergence to a poor local optimum, i.e., not discrete). The upper bounds for constraints based on equations 35 and 36 are set based on a reference value obtained by applying a small (≈ 1%) perturbation to the elasticity modulus of solid material. For clarity, it is worth considering a different view of the fully coupled step, which could be stated as finding surface displacements u * that minimize the objective function, subject to the optimization constraints, and to the additional one that ∃ρ : u * = S(ρ, λ * ). In other words, the outputs of the fully coupled step are (feasible) displacements and tractions (not the material distribution) which then become inputs of the inverse design step, tasked with producing a completely new material distribution that is discrete. An iterative procedure with the two steps could be proposed, in which case explicit penalization of non-discreteness would likely be required to maintain the discreteness achieved by the inverse design step. We emphasize that at no point do we decouple the multidisciplinary and multi-point nature of the problem, in a way in the inverse design step the FSI coupling conditions become optimization constraints, whereas in the fully coupled step they are inherently satisfied.
The optimization problem can alternatively be posed in a way that avoids the discreteness and stability issues. For example, by converting it to a compliance minimization problem with aerodynamic constraints while prescribing the topology of certain key areas, like the trailing edge. Here, we deliberately keep the design space as unconstrained as possible, and deal with the aforementioned issues in a non intrusive way to fully assess the capability of topology optimization to improve aerodynamic performance across multiple operating points.

Results
We consider two numerical examples, first a benchmark problem to demonstrate that the optimization approach of Sections 3.1 and 3.2 is adequate for geometrically nonlinear solid mechanics problems. Then, we consider a passive load adaptation FSI problem that motivated the strategy presented in Section 3.3.

Verification of the topology optimization implementation
We verify the optimization methodology chosen, and the implementation of the SIMP scheme, by reproducing published results for the classic tip-loaded cantilever problem, shown in Fig. 1. We take as reference the nonlinear topology optimization results of Buhl et al. (2000), where a material with elasticity modulus of 3 GPa and the Poisson ratio of 0.4 was considered. For this test, we use the morphology-based filters proposed by Sigmund (2007) to obtain solid-void topologies, namely the close (erode • dilate) and open (dilate • erode) strategies. The objective is to minimize end compliance (W tip ) subject to an equivalent mass constraint of 0.5, that is where δ y is the vertical displacement of the point where the load is applied. As geometric non-linearities are considered,  Referring to algorithm 1, the objective is normalized by initial value, b 0 1 = 8 and r = 2, the constraint tolerance is 0.005. The value sequence for the filter parameter β ((33)) is {0.01, 1, 4, 16, 64, 200}. The initial small value makes both morphology filters equivalent to the simpler conical filter. The loose convergence criteria for L-BFGS-B are 40 iterations or a variation of objective function value less than 10 −5 . For the tight criteria, the latter is reduced to 10 −7 and the number of iterations unlimited. The domain is discretized with 10,000 square isoparametric elements, the filter radius (defining N (i) in (33)) is twice the element size but we note that the close and open filters are applied in two stages which effectively doubles the radius of the neighborhood. We have not found ramping the SIMP exponent to be advantageous for these simple cases; a constant value of 3 is used. As the loads are large, it is not practical to start the optimization with a uniform density distribution that also respects the mass constraint; therefore, we use the optimal density distribution for SIMP exponent of 1 and linear elasticity (which is obtained with the same process described above but without ever increasing the filter parameter). Figures 2 and 3 show the topologies obtained for P=60 kN with the open filter and for P=240 kN with the close filter respectively. Figure 4 shows the convergence history for both cases, and both the obtained topologies and compliance values (W tip (60 kN) = 4.36 kJ and W tip (240 kN) = 67.0 kJ) compare favorably with the reference results (4.65 kJ and 66.5 kJ respectively). The optimization process requires on average 350 function evaluations, which again compares well with other sources (e.g., Sigmund 2007), and the obtained topologies after filtering are nearly perfectly discrete.

Baseline compliant airfoil from aerostructural shape optimization
The second example consists of a flexible airfoil operating at two distinct free-stream Mach numbers, 0.25 and 0.5, but  at the same angle of attack (AOA) of 2 • with respect to the undeformed shape. Henceforth, quantities at the lower or higher speed will be super-scripted l or h, respectively. At low speed, C l l = 0.5 must be generated and the deformation of the airfoil kept below a limit; at high speed, we wish to minimize drag.
The ideal configuration for minimum drag is a symmetric airfoil operating at zero AOA. As the airfoil deforms passively and the AOA is fixed, this configuration could only be achieved by a structure with no stiffness, which would then be unable to produce the required lift at low speed. Therefore, constraining the deformation at low speed leads to higher drag at high speed, unless the internal structure of the airfoil responds in a nonlinear way. As we intend to exploit this nonlinear structural behavior to improve passive load alleviation, a baseline airfoil for which the low speed deformation limit significantly affects performance was first designed by shape optimization with different allowed values of trailing edge displacement (y max ). For a review of the FSI shape optimization capabilities of SU2, and verification of the relevant derivatives, see Venkatesan-Crome et al. (2019). The starting point for the optimization is a NACA0012 profile with 0.5 m chord, parameterized through the free-form-deformation box (FFD) shown in Fig. 5 of which 17 points are allowed to move in the vertical direction, the bottom left point is fixed to avoid translation of the airfoil. A further constraint is added to enforce that the final area be greater or equal than the initial, that is  where α are the control point vertical displacements. For the RANS simulations, the fluid grid is an O-Grid with 77,924 nodes and sufficient wall cell size for y + ≈ 1, and the radius of the circular farfield boundary is 30 chords. The fluid is considered to be ideal and standard-sea-level properties are used for the free-stream state. The solid domain is discretized with 76,800 4 node quadrilateral elements resulting in 77,875 nodes (this level of refinement is to ensure sufficient resolution for topology optimization), the inside of the hollow region close to the leading edge is clamped, the vertical section of this region is located at 5% chord. The elasticity modulus considered is 50 MPa and the Poisson ratio 0.35. Again referring to algorithm 1, a 0 1 = b 0 1 = b 0 2 = 8 and r = √ 2, the constraint tolerance used was 0.01. Convergence criteria for the optimizer are as described for the benchmark topology problem. Table 1 shows the optimized drag and lift for different values of the deformation limit, Figs. 6 and 7 show the undeformed and deformed at Mach 0.5 (in red) airfoils for the 10 mm (0.02c) and 6 mm (0.012c) deformation limits respectively.
The optimum drag increases 3.7% by reducing the deformation limit from 10 to 6 mm. We note that the deformation constraint is not met entirely by an increase in structural stiffness (which could be achieved by increasing area) but mostly by the reduction in pitching moment that results from the reflexed camber line (compare Figs. 6 and 7).

Fully coupled optimization step
We consider the 3.7% increase in drag a significant-enough trade-off between stiffness and performance. Consequently, the airfoil optimized for y max of 6 mm at Mach 0.25 in the previous investigation is taken as the starting point for topology optimization. The trailing edge displacement constraint previously used is replaced by a compliance constraint with upper bound equal to the compliance of the  Optimization results for 6 mm constraint, undeformed (black) and deformed at Mach 0.5 (red) configurations initial design. The change of constraint function is necessary since, due to the trailing edge region not being forced to be solid, it could be possible for a much more flexible airfoil to respect the local constraint by employing large amounts of camber near the trailing edge in the deformed configuration, thereby producing the required lift at the expense of increased drag. With the area constraint no longer required, the fully coupled step is stated as where W l ref is the compliance (W = u · λd ) of the shape optimized structure at Mach 0.25 (1.141 J), and ρ the design densities (before filtering) introduced in Section 3.1. The elasticity modulus is 100 MPa, double the value used in the baseline shape optimization to allow material to be removed while maintaining compliance. The outermost 3 layers of elements (3.6% of the local thickness) are prescribed to be solid, elsewhere the initial density was 0.8, the density filter radius is 2 mm (approximately 3 times the largest element size). For this case, the SIMP exponent was gradually increased, taking the values 2, 2.5, and 3 following the strategy of algorithm 1. Figure 8 shows the history (inner iterations of L-BFGS-B) of drag and lift coefficients, compliance, and penalized objective function (f ), the large spikes correspond to the increases of SIMP exponent. Figure 9 shows the final topology and Mach number contours at the Mach 0.5 condition. As expected, the material distribution is not discrete, until around 30% chord the topology is mostly solid, around the center section the topology confers little bending stiffness acting mostly  as support for the wet surfaces. Notably the trailing edge region has compliant hinges which allow it to work as a mechanism, and the importance of this will be explained below.
The drag coefficient is reduced to 0.008812 (3.1% reduction) which is only 0.47% higher than what was obtained for the shape optimization with y max of 10 mm, but as expected the topology obtained in this step is not discrete.

Inverse design step
The inverse design step is stated as a mass minimization problem with constraints on the target deformed shape of the airfoil (35) and on the stability metric (36), that is The results for this step are shown in Fig. 10, where the black line shows the target deformed surface, and the red line shows the verified deformed surface obtained from the coupled FSI analysis loop for the topology obtained in this step. To reduce the computational cost of the inverse design step, the stability constraints based on (36) were only enabled after the design became feasible with respect to the geometric constraints ( ). Table 2 summarizes the twostep topology optimization process listing the aerodynamic coefficients and the equivalent mass (for SIMP exponent of 3) at the major checkpoints. Fig. 10 Results of inverse design step, density contours, target (black) and actual (red) deformed surfaces With the inverse design step, some of the improvement obtained in the fully coupled step is lost as the drag coefficient increases 0.39%; however, the equivalent mass is reduced by 37%. The reduced performance is mostly due to strong coupling effects that amplify the effect of the small discrepancy between target and obtained shapes. In the inverse design step, the surface error metric is constrained to 1%; however, when the performance of the resulting topology is verified, this metric increases to 2% at Mach 0.5. One of the responsible physical mechanisms is that regions of unsupported airfoil skin tend to form bumps as the airfoil flexes, and these bumps cause a reduction of the pressure (due to a local acceleration of the flow) which in turn increases the bump size. This is the main effect the stability function helps to mitigate; Fig. 11 shows how this function leads to material being added to the skin to increase its bending stiffness. Coupling a bump-based shape parameterization method (e.g., Hicks-Henne) with the topology variables could potentially reduce the importance of this positive feedback mechanism. That was not deemed necessary in this work due to the subsonic speeds; however, the effect would be more significant at transonic speeds as bumps can produce shocks.
The seemingly insignificant increase in surface error metric (from 1 to 2%) results in a 9% increase in lift (but no significant changes to the flow field). We hypothesize that this sensitivity of the design is not due entirely to our proposed method but also due to the performance metric we sought to optimize, which results in a system with strong nonlinear characteristics as its apparent stiffness changes significantly between low and high speeds. Finally we note that topology optimization results often require some form of post-processing to which the results may likewise be sensitive (for example, to remove vestigial

Comparison with conventional approaches to encourage solid-void topologies
Two common ways to encourage solid-void solutions, in problems that do not necessarily benefit from them, are to penalize intermediate densities more severely (e.g., using higher values of SIMP exponent) and to manipulate the problem formulation such that overall stiffness becomes important (e.g., embedding mass reduction into the optimization goals). The results of both these approaches are presented and discussed in this section.

Weighted objective function
First note that the low speed lift and compliance constraints require the airfoil to have a minimum stiffness; recall that the AOA is set relative to the undeformed configuration and so an airfoil that is too flexible will not produce the target lift. Therefore, with a weighted objective function of mass and drag, a stiffness-based problem can be recovered by giving no importance to drag and focusing solely on mass Fig. 13 Weighted average approach results, filtered structural density and Mach number contours, deformed configuration from fully coupled step in red for reference minimization (which would not be desirable). To test the weighted objective approach, we used weights of 0.8 for drag and 0.2 for mass (after scaling the functions), which based on the drag-mass trade-off from the fully coupled step to the inverse design step (see Table 2) should be sufficient. Moreover, numerical experiments with higher weighing of mass were less successful due to poor stability of the optimization.
The optimization process is as described for the fully coupled step except for the change of objective function and the SIMP exponent which was not ramped, being fixed at 3 from the start instead (ramping it made little difference in the results of Section 4). The convergence history is shown in Fig. 12, the optimization stalls relatively early (the last 8 iterations required 66 function and gradient evaluations) resulting in the topology of Fig. 13, where the deformed configuration from the fully coupled step is also shown for comparison.
The response of the structure is similar and the drag is lower but within one count of what was obtained with the two-step approach. The equivalent mass is 0.485, lower than in the fully coupled step, as expected, but higher than in the inverse design step, due to the less discrete material distribution (see Fig. 13).

Higher SIMP penalization
It is known that challenging topology optimization applications may require a SIMP exponent higher than that suggested by the Hashin-Shtrikman bounds to converge to solid-void solutions. To test this approach, the weighted objective optimization was continued with the SIMP exponent increased from 3 to 4 in steps of 0.25 every outer iteration of the exterior penalty method. The resulting material distribution, shown in Fig. 14, is almost indistinguishable with no significant topological changes. The drag coefficient was not improved and the equivalent mass increased to 0.542 as a result of the stiffness reduction in intermediate density areas. The lack of change could be due to the solution for SIMP exponent of 3 being a local optimum. However, the general features of the solution, dense truss-like structure at mid chord and large voids towards the trailing edge, develop very early in the optimization (note the quick decrease in mass in Fig. 12). As material is removed mostly from low strain energy areas, we hypothesize these solution features are inherent to the presence of the mass objective, and how quickly and easily it can be targeted by the optimizer (the mass function is linear).

Gradual mass minimization
One would then expect that gradually introducing the mass objective could avoid rapid convergence to the poor local optima described above, this was attempted for the results of the fully coupled step in two steps. First by adding a constraint to the optimization problem (39) to gradually lower the equivalent mass to 0.5, the initial penalty function parameters were selected such that the initial penalization was equivalent to two drag counts. Then by switching to constraining drag below 0.00885 while minimizing mass, this is required since once both mass and compliance constraints are satisfied there is no longer an incentive to improve discreteness (one of the two must be a scarce resource). Figure 15 shows the convergence history,  illustrating the more gradual decrease in mass, and Fig. 16 the obtained topology.
Although mass is reduced without significantly increasing drag, and both values are lower than those obtained after the inverse design step, the resulting topology is still far from discrete. We note that the optimization does not converge fully as the gradient of the penalized objective function (see Fig. 17) is not zero. Nevertheless, the line searches fail due to the much larger gradients at the weakly connected regions around 80% chord, thus resulting in a poor search direction.

Discussion
The trailing edge region is not highly stressed as bending moments are low; therefore, from a purely structural perspective, it does not require large regions of solid material. However, from the aerodynamic standpoint, the trailing edge is responsible for most of the load alleviation. While the topology features close to the leading edge mostly confer stiffness to the airfoil, the trailing edge behaves like a mechanism, one that under passive actuation notably increases the camber near the trailing edge at the higher speed. This localized camber leads to a more aft-loaded pressure distribution (see Fig. 18 light color line) that counteracts the effect of the reflexed camberline of the undeformed airfoil.
It is therefore plausible that this interference between aerodynamic and structural objectives leads to an intermediate design and search direction that stall the optimization. Although the problem is posed such that it benefits globally from high stiffness to mass ratio, locally (near the trailing Fig. 17 Gradual mass minimization approach, contours of derivative of penalized objective function with respect to density Fig. 18 Pressure coefficient distribution at Mach 0.5 for the baseline airfoil (6 mm constraint) and for the fully coupled step result edge) that is not what minimum drag requires. Different strategies can potentially be used to mitigate the interference between objectives without decoupling them; we note however that a strategy akin to the gradual mass minimization approach above will nearly double the computational cost, whereas the inverse design only adds 15% to the cost of the fully coupled step.

Concluding remarks
We have demonstrated how density-based topology optimization can be used to design the internal structure of a compliant airfoil with the objective of improving load alleviation. Better results were obtained than with shape optimization alone, as the nonlinear structural behavior introduced by the topology allowed the airfoil to better adapt to the different fluid loads at different speeds.
The proposed two-step methodology has addressed some shortcomings found when applying density-based topology optimization to non stiffness-based designs. In the first fully coupled step, coupled FSI simulations are considered but mass is not explicitly included as an objective or constraint. We have found this to improve the convergence of the optimization with only minor impact in performance, as it prevents the optimizer from making rapid adjustments that can cause the structure to buckle, and avoids early stopping due to convergence to poor local minima or stall due to poor search directions.
As mass is not considered in the fully coupled step, the resulting topology will, in general, not be discrete. Therefore, we have introduced a second inverse design step where a discrete topology is sought for which the airfoil response is the same. The inverse design step is computationally cheaper as it does not rely on coupled FSI simulations, and to avoid critical structures, we used a stability metric (36) whose computation requires only one FSI step (for which good initial conditions are available). Obtaining completely discrete topologies requires elaborate filtering strategies; these complicate the convergence of the optimization as the total number of iterations increases due to the need to ramp filter parameters. Moreover as the best filtering strategy is application dependent (even for simple problems), the optimization may have to be repeated for different settings. The two-step process greatly reduces the computational cost of testing different density filters at the expense of some performance, since in general only a very refined discrete topology could perfectly replicate the response of the optimum structure obtained in the fully coupled step. An iterative process, alternating between both steps, can also be considered for more complex problems but was not found necessary here. The methodology was compared with three more conventional approaches of encouraging solid-void topologies that all failed to provide a discrete material distribution for the example load adaptation problem.
We observed that due to the alleviation objective, the aerodynamic performance of the structure obtained in the inverse design step is sensitive to slight perturbations of the external shape. Therefore, it is likely that the performance would also be sensitive to approximations such as converting the not perfectly solid-void topology (due to the filter properties) to one that is completely discrete, and sensitive to eventual manufacturing inaccuracies. While repeating the inverse design step using different filter methodologies is not computationally expensive, thorough analysis of the robustness of the design with respect to any inaccuracies would still be required. Finally as nonlinear modeling of the structure is considered and strains are large, load-path analysis would have to be conducted for the complete operating range.
Overall the proposed approach offers a route for RANSbased topology optimization of FSI systems with focus on aerodynamic performance rather than on the realizability of the resulting topologies.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.

Replication of results
The source code developed for this work, the data, and the scripts required to run the optimizations are available at https://github.com/pcarruscag/SU2 under the tags SMO full coupled and SMO inv design.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.