Topology Optimization of Fluid-Structure Interaction Problems with Total Stress Equilibrium

This work extends force coupling in the topology optimization of fluid-structure interaction problems from hydrostatic to total stresses through the inclusion of viscous stress components. The superconvergent patch recovery technique is implemented to remove the discontinuities in velocity derivatives over the finite elements boundaries. The sensitivity analysis is derived analytically for the superconvergent patch recovery approach and further verified through the use of the complex-step derivative approximation method. Numerical examples demonstrate a differentiation in the optimized designs using pressure vs. total stress coupling depending on the flow characteristics of the design problem.


Introduction
Fluid-structure interactions (FSI) refer to problems where a deformable or a movable structure interacts with a flowing or stationary fluid.Examples of FSI in everyday life are abundant; the flow of blood through flexible arteries, the deflection of an elastic airplane wing under aerodynamic loads, and the vibration of bridges and tall buildings under wind loads [1].A large body of research work has been dedicated to the numerical analysis of these interactions [2,3], and most multiphysics simulation packages -commercial and open source -include some degree of capability for solving these problems.However, the application of numerical design tools such as topology optimization (TO) to these problems is still lagging behind.This is mainly due to the multidisciplinary nature of the combined topology optimization of f luid-structure interaction (TOFSI) problem and its strong nonlinear behavior resulting in numerous stability and convergence issues.
Density-based TO methods first appeared in the seminal work by Bendsoe [4] which marked the start of modern TO in its current form.In density-based methods, each discretization unit (e.g.finite element) is interpolated between material and void throughout the optimization process.Penalization of a property of interest is used to force the final solution towards a discrete 0/1 state [5].Later, TO was extended to Stokes fluid flow in [6] using a parameter to control permeability for design parametrization.In [7], the authors extended the work to the Navier-Stokes equations and used the analogy of a 2D channel flow with varying thickness for design parametrization.They also recognized, along with Evgrafov [8], the similarity between this design parametrization and Brinkman equations [9] of fluid flow in porous media [10, p. 15].This approach -later termed Brinkman penalization -became the de facto method for TO of fluid flow without much regard to its physical interpretation.
The first application of TO to FSI problems appeared in [11].The authors used a segregated domain formulation to define the three-field FSI problem (i.e.fluid, solid, and mesh).In the segregated formulation, there is no overlap between the fluid and the solid computational domains thus only allowing for dry TO (i.e. the fluid-structure interface is not optimized).This formulation was used in [12] and extended to mechanism design in adaptive wings in [13].In [14,15], the authors utilized the three-field, segregated domain formulation in the design of wing skeletons of micro air vehicles using a two-material formulation; wing skin or carbon fiber reinforcement.Although -technically speaking -this was a form of wet TO, due to using a two-material formulation there were no fluid porous elements, and hence from a force coupling perspective it's similar to the dry TO described in [11].In [16], the authors performed TO on a flexible micro-fluidic device, however, they only considered the effect of structural deformations on the fluid flow and ignored the fluid forces on the structure.In [17], a polytree-based adaptive polygonal finite element method was used for TO of fluid-submerged breakwater barriers.However, they used a low fidelity model such that the initial fluid loading was used to calculate the forces without further updates based on the changes in the design domain.Hence the fluid forces were calculated on the fully solid structure and no porosity was incorporated in the fluid domain.
The first appearance of the unified domain formulation was in [18] where the fluid and solid computational domains overlapped such that the whole design domain -including the fluid-structure interface -could be optimized, hence the term wet TO.In [18], the author extended the fluid-solid force coupling to the entire fluid domain including porous elements but only considered the hydrostatic component of the fluid stress tensor.The divergence theorem was utilized to transfer the integral from surface to volume.The same unified formulation and force coupling was later utilized in [19,20,21] by the same author.More recently, in [22], the authors revisited the density-based unified formulation with special emphasis on the multidisciplinary interaction between the fluid and the solid and its effects on the optimized designs.They also investigated different fluid flow conditions and demonstrated the significance of careful selection of interpolation functions and parameters to ensure convergence.Nonetheless, all these works were limited to hydrostatic traction equilibrium.
In the present work, we aim to extend the traction equilibrium condition to include the total stress tensor such that viscous stresses are considered in addition to the previously-considered pressure stresses.The remainder of the paper is organized as follows; in section 2, we state the governing equations for the fluid-structure interaction problem.In section 3, we detail the finite element formulations of the governing equations and the traction equilibrium condition as well as the relevant sensitivity analysis.In section 4, we present the test problem to be used in the following sections.In section 5, we discuss the implications of the discontinuities in velocity derivatives at the elemental boundaries and how to overcome this issue.In section 6, we demonstrate with numerical examples the effect of using pressure vs. total stress force coupling.Finally, in section 7, we list the conclusions.

Governing Equations of the Fluid-Structure Interaction Problem
As further detailed in Sec. 3, we limit our discussion to infinitesimal strains in this work.For large deformations, a mesh deformation technique needs to be implemented and the fluid governing equations need to be redefined in the deformed configuration, cf.[18, p. 594] for such a derivation in a TOFSI context.In our implementation of the unified domain formulation, the fluid computational domain Ω f spans the whole computational domain Ω .The solid computational domain Ω s is fully contained within the fluid computational domain.The solid non-design domain Ω n , if any, along with the design domain constitute the solid domain.Such that: In Ω d , a design variable ρ per each finite element determines whether it's pure fluid, pure solid, or in between.More specifically, ρ represents the solid density of the design element; i.e. ρ = 1 for pure solid, ρ = 0 for pure fluid, and 0 < ρ < 1 for porous solid.The solid/fluid interpolation required in density-based TO is achieved by: (i) appending the fluid momentum equation with a Brinkman penalization term −α(ρ)v as a volume force where the inverse permeability1 α(ρ) is dependent on the design variables ρ, and (ii) interpolating the structural stiffness E(ρ), i.e.Young's Modulus, as a function of the design variables ρ.
The Brinkman penalization term enforces zero fluid velocity for ρ = 1, porous flow conditions for 0 < ρ < 1, and pure fluid flow for ρ = 0.This way solid and fluid are allowed to co-exist within the same finite element, which is necessary to transform the discrete optimization problem into a continuous one.If a proper physical interpretation of the Brinkman penalization concept is sought, we assume that the porous media consists of a volumetric mixture of fluid and solid phases, the solid phase is homogeneous and isotropic on the macroscopic scale, and all the pores are interconnected [10, p. 4].
For the fluid flow, the Navier-Stokes equations are used in their viscous, incompressible, steady-state form.The strong form of the partial differential equations is as follows [23, p. 10]: where v is the fluid velocity, ρ f is the fluid density, σ f is the Cauchy fluid stress tensor, f f is the external fluid force, p is the hydrostatic pressure, and µ is the fluid dynamic viscosity.In Eq. 7, the inverse permeability α(ρ) is dependent on the Brinkman penalization upper and lower limits, α max and α min respectively, and the Brinkman penalization interpolation parameter p α that is dependent on the problem physics; e.g.Reynolds number [22, p. 974].
As for the structure, the Navier-Cauchy equations are used assuming linear elasticity with infinitesimal deformations under steady-state conditions.The strong form of the partial differential equations in the stress-divergence form is as follows [24, p. 168]: (10) (11) where σ s is the Cauchy solid stress tensor, f s is the external solid force, C s is the solid elasticity tensor, is the infinitesimal strain tensor, and u is the solid displacement.The elastic modulus E(ρ) is interpolated using the modified SIMP approach where E min and E max are the lower and upper limits on E representing the void and bulk elastic moduli, respectively.p E is the elastic modulus penalization parameter.The essential boundary conditions are defined as follows: Fluid Outlet: Structural Prescribed Displacement: Note that the fluid no-slip boundary condition in Eq. 12 is only defined on the external domain boundaries aside from the inlet and outlet such that ∂Ω f = Γ v0 ∪ Γ vin ∪ Γ vout .The volume force term appended to the fluid momentum in Eq. 5 automatically enforces a no-slip condition wherever needed within the solid domain and/or its fluid-structure interface.The natural boundary condition in an FSI problem is defined as follows: where n f and n s are the normals to the fluid and solid surfaces, respectively.After calculating the fluidic forces per each finite element and before the global assembly, the forces must be multiplied by the pressurecoupling filtering function defined as [18, p. 607]: where Υ max and Υ min are assumed 1 and 0 respectively.p Υ is the pressure-coupling filter function interpolation parameter.The main intention of using the pressure-coupling filtering function is to prevent severe deformation of low density elements and potential singularities in solving the structural finite element system.While the fluid-structure interface Γ FSI is explicitly defined in pure fluid-structure interactions (i.e.no porous media), this is not the case in density-based TO problems where all solids are porous to some degree.Hence, to apply the traction coupling condition in this formulation properly, Γ FSI must be taken as all edges of all solid finite elements.More discussion on this condition and how to apply it in a finite element form follows in subsection 3.2.
In the next section, we transform the strong form of the governing equations into the discretized finite element form and extend the traction equilibrium condition described generally in Eq. 16 to include viscous stresses.
3 Finite Element Formulations

Fluid and Solid Governing Equations
In this work, conformal or matching discretization is assumed at the fluid-structure interface to ensure simple, yet accurate, force coupling.For non-conformal meshing, other mapping techniques can be used to enforce momentum conservation at the interface [25].The discretized finite element form of Navier-Stokes equations are obtained by multiplying the strong form of the continuity and momentum equations (Eqs.4 and 5) with suitable weight functions and integrating the resulting expressions over the appropriate computational domain.The divergence theorem (i.e.integration by parts) is typically implemented to reduce the continuity requirements on the interpolation/shape functions used to approximate the solution by moving the differentiation to the weight functions [23].
In this work, we implement the standard Galerkin method of weighted residuals where the interpolation/shape functions themselves are used as test/weight functions.The resulting finite element model is of the velocity-pressure (or mixed) type where both velocity and pressure are unknowns and are solved for simultaneously.We utilize a quadrilateral meshing which is regular and structured wherever possible.To satisfy the Ladyzhenskaya-Babuska-Brezzi condition (cf.[23, p. 176]), Q2Q1 Lagrangian finite elements (i.e. 9 velocity nodes and 4 pressure nodes) are used with a low to moderate Reynolds number to avoid using stabilization techniques that artificially -but not necessarily accurately -dampen the discontinuities.The continuity equation (Eq.4) is typically weighted by the pressure shape function Φ while the momentum equation (Eq.5) is typically weighted by the velocity shape function Ψ, both are in column vector form.The resulting finite element system in 2D on the elemental level is as follows: The coefficient matrices in the finite element form are defined as ( ´+1 −1

´+1
−1 and dξ dη are implied)2 : where v1 , v2 , and p are the nodal velocities in x and y and nodal pressures, respectively.|J| is the Jacobian determinant and ξ and η are the natural coordinates.F f 1 and F f 2 are externally applied nodal fluid forces which are not used in this work since the fluidic boundary conditions (Eqs.12 to 14) are implemented directly by setting nodal velocities/pressures to their appropriate values.Appropriate global assembly of Eq. 18 is implemented and the resulting nonlinear system is solved using the undamped Newton-Raphson method [23, p. 190].
As for the Navier-Cauchy equations, their discretized finite element form is obtained using the principle of virtual displacement following the procedure in [26, p. 153].To simplify the force coupling, we utilize Quadratic Lagrange finite elements (i.e. 9 displacement nodes) for the solid domain.In 2D, the resulting system on the elemental level is as follows: where û is the nodal displacements and F s is the nodal solid forces.The solid stiffness matrix k s for plane strain conditions is obtained following the procedure in [26, p. 164] with the only difference being the dependence of the elastic modulus on the design variable as detailed in Eq. 11, and appropriate global assembly is also performed to construct the global solid stiffness matrix.Isoparametric elements are used for both fluid and solid domains, so the geometry is described by the same shape functions describing the velocities/displacements.

Traction Coupling
In literature, the traction equilibrium condition was derived considering hydrostatic stresses only (cf.[18, p. 601] and [22, p. 989]), which is suitable for some cases with negligible viscous stresses.In this work, we extend the traction equilibrium condition to include the viscous stresses as well.The fluidic forces applied on the solid structure are calculated directly from Eq. 16 as follows: which can be put into the finite element form through multiplying by the velocity shape -and test -functions vector Ψ then discretizing: After transformation from global to natural coordinates, Eq. 25 can be further detailed as follows: where the normal vector components n f x and n f y can be calculated from elemental nodal coordinates, and dξ s designates the natural coordinate along the finite element edge of interest in 2D.Note that J s is the surface Jacobian that is different from the typical area Jacobian J (cf. [23, p. 97] or [26, p. 356]).
In previous works on TOFSI problems (and in works on FSI with porous media in general), the divergence theorem is typically used to transform the surface integral into a volume integral as follows (cf.[18, p. 598] and [22, p. 989]): By applying the divergence theorem in Eq. 28 to Eq. 25, the surface integral form in Eqs. 26 and 27 is transformed into a volume integral form as follows: Note the appearance of velocity first derivatives in the force coupling equations using surface integrals (i.e.Eqs. 26 and 27).In the volume integral form of the force coupling (i.e.Eqs.29 and 30), velocity second derivatives and pressure first derivatives appear due to the use of the divergence theorem.The appearance of these derivatives and their effect on traction equilibrium when using total stress definitions is the topic of discussion in section 5.

Sensitivity Analysis
Typically, in TO problems, the number of design variables greatly exceeds the number of constraints.Hence, it is better from a computational perspective to use the adjoint method in calculating the sensitivities as opposed to the direct method (cf.[27] for a discussion on the adjoint method in an FSI context).Given an objective function (or a constraint) f , the sensitivity of f w.r.t. the design variable ρ i is calculated as follows: where r is a vector of all state variables (i.e.solid displacements and fluid velocities and pressures), R represents the set of all (i.e.solid and fluid) finite element equations, and λ is the adjoint vector which is distinctive for each objective function or constraint.
It can be seen that the derivatives of the force coupling equations w.r.t. the fluid state variables are needed as a part of the ∂R/∂r term.For the case of hydrostatic coupling implemented earlier in [18] and [22], only the derivatives w.r.t. the nodal pressures are needed.These derivatives are evaluated for the surface integral form directly from Eqs. 26 and 27 as follows: and for the volume integral form, they are evaluated directly from Eqs. 29 and 30 as follows: 4 Description of Test Problems Before we discuss the details of the total stress equilibrium equations in section 5, we first describe the setup of the numerical example to be used in the following sections.
The original column in a channel problem was first discussed in a TOFSI context in [18, p. 610] and has been used later as a benchmark problem in a number of works on TOFSI.In [22], Lundgaard et al. developed a modified version 3 , where they (i) increased the relative size of the design domain with respect to the whole computational domain, (ii) enlarged the geometry from the micro to the macro scale, and (iii) generally strengthened the fluid-structure dependency.We utilize the modified version of this problem as the numerical example of choice in this work.In the remainder of this work, we drop the designation "modified" and simply refer to this problem as the "column in a channel" test problem.
As detailed in Fig. 1, the problem features a 1.4 × 0.8 m rectangular design space (light gray) placed inside a 2 × 1 m rectangular channel.To avoid trivial solutions, a 0.5 × 0.05 m non-design, elastic, solid column (dark gray) is placed within the design space to force the optimizer into a more sophisticated solution than a simple bump at the bottom of the channel.Note that even though the non-design column is not included in the design domain, it is still governed by the structural equations and is included in the calculation of the objective function.
The top and bottom surfaces of the channel Γ v0 have a no-slip condition applied.A parabolic laminar flow profile is applied at the inlet Γ vin on the left with a maximum velocity of 1 m/s, and a zero pressure condition is applied at the outlet Γ vout on the right.The bottom surface of the design and non-design domains Γ u0 is fixed to a ground structure.The fluid density ρ f is assumed 1 kg/m 3 .The maximum and minimum Brinkman penalization limits are taken as α max = 1e + 9 and α min = 0 as suggested in [22, p. 974] and further discussed in [28].The solid elastic moduli parameters E max and E min are assumed 1e + 5 Pa and 1e − 5 Pa, respectively, while Poisson's ratio ν is taken as 0.3.The fluid dynamic viscosity µ is to be varied to produce different Reynolds numbers which are calculated according to Eq. 42 as follows: where L c is a characteristic length that is taken as the width of the inlet, i.e. 1 m, in agreement with [22, p. 971].A regular, structured mesh of 45,300 (i.e.302 in x and 150 in y) quadrilateral Q2Q1 finite elements is implemented. 3Lundgaard et al. [22] called this problem "the wall".

Total Stress Equilibrium in Porous Fluid-Structure Interaction Problems
It's well-known that in the Galerkin finite element formulation with typical Lagrangian shape functions, the primitive variables (i.e.velocities and pressure in fluid flow) are continuous across element boundaries, however, no such continuity is enforced on their derivatives [29].This discontinuity in derivatives poses a challenge to the calculation of relevant derived quantities that utilize these derivatives such as fluxes andpertinent to this work -traction at the fluid-structure interface.
A number of techniques have been discussed in literature to overcome this issue.Bazilevs and Hughes [30] enforced the Dirichlet no-slip boundary condition weakly in the context of turbulent fluid flow instead of the typical strong nodal imposition.Their mean velocity profile results showed good agreement with comparable direct numerical simulation results.Oshima et al. [31] developed a consistent method as a post-processing step where the fluxes or stresses are treated as unknowns which are solved for weakly after the typical Galerkin finite element solution is obtained.Zienkiewicz and Zhu [32,33] developed a superconvergent recovery technique for calculating the derivatives of finite element solutions at nodes.In this work, we implement the discrete version of this technique to overcome these discontinuities in the calculation of total stress coupling enforced at the fluid-structure interface.In particular, we apply the technique to recover the velocity derivatives w.r.t.coordinates x and y.

Recovery of Discontinuous Velocity Derivatives
The superconvergent recovery technique by Zienkiewicz and Zhu [32,33] has its basis in the discovery of superconvergent points within each finite element where derivatives are known to have extraordinary accuracy compared to the rest of the finite element (cf.[32,33] and references therein).These points correspond to reduced integration, so for the quadratic quadrilateral finite elements used in this work (3 × 3 velocity nodes), these superconvergent points are the 2 × 2 Gaussian integration points with equal weights.The procedure we implement is described in detail in [32, p. 1335] and in a fluid mechanics context in [29, p. 952].Once the recovered velocity derivatives at the nodes are calculated, they can be interpolated within each finite element using the same shape functions used to interpolate the primitive variables, i.e. velocities.Hence, this technique lends itself well to the volume integral form of force coupling in Eqs.29 and 30 as the second derivatives of velocities are evaluated at the 3 × 3 Gaussian points within each finite element using the appropriate derivatives of the shape functions multiplied by the recovered nodal velocity derivatives as detailed later in Eqs.52 and 53.
We reiterate the recovery procedure here as it's pertinent to the sensitivity analysis discussion in subsection 5.2.We designate the recovered values with a bar on top, so the recovered nodal velocity derivatives are v1,x , v1,y , v2,x , and v2,y .The scalar recovered velocity derivative within a finite element v i,j is related to the corresponding recovered nodal values of that particular element by interpolation using velocity shape functions as in: Each nodal value of the vector vi,j is calculated within one or more 4 recovery patches where it's represented by a polynomial expansion of the following form: where P contains the appropriate polynomial terms and is a function in coordinates x and y within the patch 5 , and a i,j -a vector of unknown parameters of the same size as P -is solved for within each patch through a least square fit minimization problem as follows: and v i,j is evaluated at the superconvergent points (i.e. 2 × 2 Gaussian points for a Q2Q1 element) within each finite element in the recovery patch using the appropriate relation as in: The vector P is designated with a subscript O when evaluated at the to-be-recovered nodes as in Eq. 44 and with a subscript × when evaluated at the superconvergent Gaussian points as in Eq. 45.This distinction is critical in the sensitivity analysis discussion in subsection 5.2.
Note that within each patch, the four different derivatives -i.e.v1,x , v1,y , v2,x , and v2,y -are recovered using the same P, however, a different a i,j is calculated for each derivative, hence the subscript designation i,j for a i,j .A few remarks about our implementation of the superconvergent patch recovery technique are in order: a.Assuming P = {1, x, y, xy, x 2 , y 2 , xy 2 , x 2 y, x 2 y 2 } T for Q2Q1 elements and depending on the geometric scale of the problem, singularities might arise in solving the linear system a = A −1 b (where A = P P T and b = P v i,j ) due to how the matrix A is constructed.For instance at point (x, y) = (1, 1), all terms of P are equal to 1.To overcome this issue, we elected to scale all x and y coordinates within each patch.We arrived at the following scaling relations using -may we say -'educated' trial and error: where x and ŷ are the nodal coordinates of a random element within each recovery patch, usually the first element to be considered in the programming for loop within each patch.It's possible that these scaling relations might need to be tuned based on the geometry of each problem as this was not rigorously tested for different problems.
b. Depending on the regularity of the mesh and the location of the to-be-recovered nodes, some patches may contain less or more than 4 elements.In particular, for patches with less than 4 elements, the degree of the polynomial -hence the number of terms -used to construct P has to be reduced to avoid singularities.From our numerical experiments, patches containing two elements (e.g. the center patch node is at an external boundary) may use 4 polynomial terms, patches containing three elements (e.g. in the middle of an irregular mesh) may use 6 polynomial terms, and patches containing 4 or more elements may use all the 9 polynomial terms.
To illustrate the effect of using the superconvergent patch recovery technique, Fig. 2 shows the original vs. recovered ∂v 1 /∂x and ∂v 2 /∂y derivatives at a square area with coordinates x = [0.9m : 1.1 m] and y = [0.4m : 0.6 m] of the column in a channel test problem, defined in section 4. In solving the fluid analysis problem and extracting these derivatives, we utilized the default values described in section 4 in addition to the following parameters; µ = 1 Pa•s, p α = 5.25e − 6, p E = 1, and p Υ = 1.An initial guess of 0.1 is used for all elements in the design domain (i.e.ρ = 0.1 in Ω d ).It can be clearly seen that using the superconvergent patch recovery technique restores the continuity of the velocity derivatives across the elemental boundaries.However, the following phenomena are observed: a.The polynomial degree representing the velocity derivatives is raised from bilinear in the original derivatives to biquadratic in the recovered ones.
Original @u=@x -15 -10 b.Some peaks in the original velocity derivatives are smoothed out in the recovered derivatives as evident at the two top corners of the non-design column, i.e. at (x, y) = (0.975, 0.5) and (1.025, 0.5) m.
c. Due to the non-design column having almost zero velocities, the velocity derivatives within are almost vanishing and don't make a lot of sense.
Next, we move on to rederive the force coupling equations using the recovered velocity derivatives.The modified version of the surface integral form of force coupling with recovered derivatives is as follows: and the modified version of the volume integral form is as follows:

Sensitivity Analysis using the Recovered Velocity Derivatives
As detailed in the discussion on the sensitivity analysis (i.e.subsection 3.3), the derivatives of the solid governing equations w.r.t. the fluid state variables are needed as a part of the ∂R/∂r term.Before we discuss the sensitivity analysis any further, we share a thought on the discontinuities in the forces sensitivities.
Even though the derivatives of the forces w.r.t.velocities are independent of the nodal velocity values, there are still discontinuities at the elemental boundaries.By taking a closer look at any velocity derivative, take for instance ∂ ∂v1 /x = ∂ ∂Ψ /x T v1 , evaluated at an elemental boundary from its two sides, the discontinuities arise due to a combination of (i) differing velocity components in v1 6 , and (ii) differing shape functions derivatives, i.e. components of ∂ ∂Ψ /x.With the nodal velocities missing from the forces sensitivities, we focus on the discontinuities in the shape functions derivatives.The reason for discontinuities in the shape functions derivatives lies in the irregularity of the mesh and how this affects the Jacobian matrices of adjacent elements.In a highly regular mesh, these discontinuities almost disappear and the sensitivity analysis of the forces may be calculated directly by differentiating the original equations before recovering the velocity derivatives, i.e. by differentiating Eqs. 26 and 27 for the surface integral form and Eqs.29 and 30 for the volume integral form.However, in the interest of generality in this work, we derive the sensitivity analysis equations after recovering the velocity derivatives.
The sensitivity analysis of the force equations w.r.t.velocities in the surface integral form is as follows: while the volume integral form is as follows: The derivatives of a recovered nodal value vi,j w.r.t. the nodal velocities of a specific element ve i are calculated, in vector form, as follows: with no summation on i in the LHS term.m is the number of patches sharing a particular recovered node, n is the number of superconvergent points within a particular patch l, and q is the number of superconvergent Gaussian points within a particular element e which is 4 for Q2Q1 elements.Note that the use of superconvergent patch recovery extended the dependency between the adjacent elements' nodes, further complicating the sensitivity analysis.However, fortunately, Eq. 62 is independent of nodal velocity values, hence it need not be updated for new state variables.Nevertheless, in cases with mesh deformation, the equation needs to be solved after any changes in nodal x and y coordinates.

Verification of the Sensitivity Analysis
In this subsection, we utilize the complex-step derivative approximation method (cf.[34]) to verify the sensitivities calculated for the pressure as well as total stress coupling cases.The complex-step derivative where Im(f ) indicates the imaginary part of f .Utilizing the column in a channel test problem defined in section 4, we randomly sampled eight finite elements within the design domain and calculated the sensitivities at these elements using the complex-step derivative approximation method against the analytically-derived sensitivities discussed in subsection 5.2. Figure 3 shows the location of the finite elements of interest within the design domain.For the results in table 1, we utilized the parameters used in subsection 5.1 without filtering.The error in table 1 is calculated by subtracting the analytical sensitivity from the complex sensitivity then dividing over the analytical sensitivity.As can be noted, the errors are within a reasonable range attesting to the validity of the derived analytical sensitivity.

Numerical Experiments
In this section, we present some numerical results through solving the test problem described in section 4. Our main intention in the following discussion is to demonstrate the difference between the pressure vs. total stress force coupling.The method of moving asymptotes (MMA) is used to update the design variables [35].In this work, all the derivatives are derived analytically -by hand -and the entire TOFSI problem is coded and solved in MATLAB.
In general, TOFSI with pressure coupling uses less computational resources when compared with TOFSI with total stress coupling.As mentioned in the last paragraph in section 5.2, in the absence of mesh deformations, the derivatives of the fluidic forces w.r.t. the state variables are calculated only once in the first iteration.In the subsequent iterations, the pressure-coupling filtering function is updated with the new design variables and multiplied by the derivatives of the fluidic forces w.r.t. the state variables.Based on our coding practices, we estimate the increase in computational costs to be below 10%.
The column in a channel problem (described in section 4) is solved with the objective function of minimizing the compliance of the structure at a specific volume fraction as follows: where V i is the volume of finite element i, V 0 is the volume of the entire design domain, and the volume fraction V f is set to 0.1 following [22, p. 976].The volume fraction is further used as the initial guess for all the optimized designs.The problem is optimized for two Reynolds numbers; Re = 1 and 10, obtained by varying the fluid dynamic viscosity as discussed in section 4 such that µ = 1 and 0.1 Pa•s, respectively.

In Search of Discrete Designs
In [18], Yoon used constant interpolation parameters to solve the original column in a channel problem.Specifically, he used p E = 3, p Υ = 3, and p α = 11 for different Reynolds numbers, albeit the Brinkman penalization interpolation function was defined differently from Eq. 7 and more similar to the typical SIMP equation.Figure 4 shows a comparison between the Brinkman interpolation function used by Yoon [18] and the interpolation function implemented in this work (Eq.7), showing that the former is generally more sensitive to intermediate density elements.Later in [22], Lundgaard et al. solved the modified column in a channel problem using the robust formulation [36] with p E = 1, p Υ = 1, and p α calibrated for different Reynolds numbers.Lundgaard et al. [22] brought considerable attention to the effect of p α on the optimized designs.
Initially, we attempted to solve the problem using continuation on p E , p Υ , and p α in order to avoid the exacerbated computational costs associated with using the robust formulation.However, TOFSI problems proved to be notoriously non-convex and whatever computational costs were saved by avoiding the robust formulation, were consumed in the numerous trial-and-error attempts to calibrate the interpolation parameters.Moreover, the best combination we could find still contained considerable intermediate densities.Eventually, we settled on a combination of the robust formulation and continuation on p E and p Υ .From our perspective, obtaining close-to-discrete designs takes precedence over enhancing the objective function.
To generate the following designs, the problem is started with p E = 1, p Υ = 1, and β = 4.These parameters are updated such that at iterations {21, 41, 61, 81}, p E is raised by {0.5, 0.5, 1.0, 1.0}, p Υ is raised by {0.5/δ E|Υ , 0.5/δ E|Υ , 1.0/δ E|Υ , 1.0/δ E|Υ }, and β is set to {8, 16, 32, 64}.δ E|Υ is a parameter that is used to modify the penalization of p Υ against that of p E , hence it takes values greater than or equal to unity.For δ E|Υ = 1, p E and p Υ experience the same penalization.As δ E|Υ increases, p Υ gets less penalized than p E .For the robust formulation, we followed Eqs.10-13 in [22, p. 973] with projection thresholding values of η n = 0.5, η d = 0.49, and η e = 0.51 for the nominal, dilated, and eroded designs7 , respectively.The relative filtering radius is set to 5.3, which is a bit high but is necessary to avoid thin members that tend to contain intermediate density elements even at high β values.
To achieve this relatively low number of iterations per step, we used the default settings of the MMA code in its min/max form [37, p. 3] except for the following changes; (i) the move limit is set to 0.1, (ii) a positive offset of 1 is added to the objective function before inputting into the MMA code to increase the constraints violation, and (iii) the objective function (and its derivatives) is multiplied by 100 to raise the order of magnitude of the derivatives for the cases of Re = 10.These modifications generally increase the "aggressiveness" of the MMA optimizer with a slight detriment to the optimality of the final design.In our numerical experiments, the maximum change in the design variables is generally above 0.05 and typically higher at the start and after each continuation update.The volume fraction constraint is applied on the dilated design with its value updated every 2 iterations to account for this increased aggressiveness, cf.Eq. 14 in [38, p. 774].
From our numerical experiments, we didn't notice any convexification effects in implementing continuation of p α , hence it's set to a single value throughout the entire optimization process.However, we noticed that the optimized designs are far more sensitive to the selection of p α than to the selection of p E or p Υ .This is in agreement with the conclusion reached in [22, p. 990].To settle first the selection of δ E|Υ , hence p Υ , we present the optimized designs for different p Υ for two p α values for the case of Re = 1. Figure 5 8 shows the effect of different δ E|Υ on the optimized designs under pressure coupling for p α = 5.25e − 6.It seems that, generally, close-to-discrete designs are better obtained for δ E|Υ somewhere between 2 and 3. Meaning that intermediate density elements are exposed to a relatively stronger fluid forces, hence driving their density values towards the discrete limits.In addition, the effect of δ E|Υ is far more significant when p α is not properly selected.A value of p α = 2.5e − 6 seemed more appropriate for the flow conditions of Re = 1 for both pressure and total stress coupling as can be seen in Figs. 6 and 7, respectively.It can be noted that the effect of δ E|Υ in the case of the more appropriate p α = 2.5e − 6 is far less significant, yet still noticeable.Similarly to the case of p α = 5.25e − 6, a value of δ E|Υ between 2 and 3 produces better designs for the case of p α = 2.5e − 6 as well.It's critical to mention that the above discussion on the results presented in Figs. 5, 6, and 7 concerns the effect of the relative values of p Υ vs p E rather than the absolute values of p Υ .
To further elaborate on the effect of p α on the optimized designs, we solved the case of Re = 10 using both pressure and total stress coupling for a range of p α values.A value of δ E|Υ = 2 is used such that at iterations {21, 41, 61, 81}, p E is set to {1.5, 2.0, 3.0, 4.0} and p Υ is set to {1.25, 1.50, 2.00, 2.50}.Some interesting conclusions can be drawn on the effect of p α : a.One important remark we noted in this work and in a previous work [28], is that the pressure is much more sensitive to the selection of p α than the velocity.Meaning that a higher p α is needed to get an accurate pressure field than is needed to get an accurate velocity field.By "accurate" we mean in comparison to a segregated domain, body-fitted formulation approach.
b. Large p α means that the pressure is "over-sensitive" to intermediate densities.Hence, the system tends to "abuse" this sensitivity and build intermediate density "walls" that redirect the flow and diffuse the pressure around the main structure.Thus reducing the objective function as the structure experiences a non-physically lower force than it should.This can be seen in Figs.8a to 8e for the pressure coupling scenario and similarly in Figs.9a to 9i for the total stress coupling scenario.
c. Small p α means the pressure is "under-sensitive" to intermediate densities.The system then "abuses" this insensitivity to pressure and extend thin structural members further to the left in a way that     better supports the non-design column, which is still experiencing the correct fluidic forces.These thin structural members are "immune" to pressure loads due to the small p α selected thus rendering the design non-physical and wouldn't pass a thresholding test.This phenomenon can be observed in the pressure coupling scenario in Figs.8l to 8o.
At this point, it makes sense to formally define the appropriate p α for a certain case as the value that doesn't result in non-physical, intermediate density walls and at the same time resolves the correct pressure field.Thus taking pressure loading into account during the design process and potentially passing a thresholding test.Building upon the above three remarks, it seems that the appropriate p α for pressure coupling tends to be relatively higher than the appropriate p α for total stress coupling.That is why the total stress coupling cases tend to display the same behavior as the pressure coupling cases towards changing p α , albeit at a lower p α .Examples of this observation are: a.The total stress coupling cases in Fig. 7 still show intermediate density walls at δ E|Υ = 1 & 1.5 in contrast to the pressure coupling cases (Fig. 6) not showing intermediate density walls even at δ E|Υ = 1 for the same p α .
b.The same observation is more pronounced in Fig. 8 where the pressure coupling cases stop demonstrating intermediate density walls at p α = 5e − 7 while the total stress coupling cases (Fig. 9) still demonstrate intermediate density walls down till roughly p α = 6e − 8.
Based on these observations, we recommend -for calculating the appropriate p α for a certain problem -to start from a high enough p α that results in intermediate density walls, then reduce its value in steps until they disappear without any further reduction.The pressure field should be also monitored visually to ensure it can still "see" the structure at the lower p α values.

Effect of Changing the Reynolds Number on the Optimized Designs
Changing the Reynolds number seems to have a similar effect to what Lundgaard et al. [22, p. 978] described.Increasing Re tends to shift more material towards the upstream half to better strengthen the structural member on this side which faces increasing forces from the fluid flow.The same phenomenon is observed for both pressure and total stress coupling scenarios.equations.We noted that the effect of Re cannot be isolated from its constituents (i.e.viscosity, velocity, density, and length); hence how Re is changed was concluded to be more significant than the value of Re itself.It is not straightforward to extend these conclusions to areas of intermediate and low density due to the strong non-linearity of the fluid flow equations, but we believe it's worth investigating in future work.
To check this hypothesis, we solved two examples in which we used v max = 10 m/s and µ = 1 Pa•s to get Re = 10.This is in contrast to the previous Re = 10 designs (Figs. 8 and 9) in which we used v max = 1 m/s and µ = 0.1 Pa•s.We used p α = 5e − 7, δ E|Υ = 2, and 0.01 as a scaling factor for the objective function and its derivatives.Interestingly, the pressure coupling design (Fig. 10a) bears strong resemblance to the pressure coupling design obtained for p α = 5e − 8 while using µ to change Re (Fig. 8o).Similarly, the total stress coupling design (Fig. 10b) bears strong resemblance to the total stress coupling design obtained for p α = 5e − 8 while using µ to change Re (Fig. 9o).This confirms our hypothesis that Re on its own is not enough to completely characterize fluid flow behavior in FSI problems.

Effect of Pressure Spikes
In the existence of external sharp corners, pressure spikes -in the form of extreme high and low localized values -are prone to appear.For instance, Fig. 11 shows the velocity and pressure fields at the first iteration of the Re = 1 case, where these pressure spikes can be observed at the top sharp corners of the non-design column.Note that the type of force coupling used doesn't affect the resolved velocity and pressure fields, but rather how the fluidic forces are extracted after.Hence, the velocity and pressure fields of the initial design shown in Fig. 11 are exactly similar for either type of force coupling.The sensitivities, on the other hand, are indeed affected by the type of force coupling.
In the pressure coupling cases, the sensitivities are greatly affected by these pressure spikes, which led the optimizer to place material in this vicinity to alleviate their effect.However, in the total stress coupling cases, the sensitivities were not affected as much by these pressure spikes, and the available material was used somewhere else to better support the design.This is evident in Fig. 12, which shows the sensitivities for Re = 1 at the second iteration for the pressure vs. total stress coupling cases.
A pressure vs. total stress cross-check is performed in table 4 by analyzing each design for pressure as well as total stress coupling (the yellow highlight indicates the best performance).It is clear that each design performs well for its intended purpose and that excluding the viscous stresses may alter the optimized result depending on the flow characteristics.

Conclusions
In this work, we extended the force coupling in the topology optimization of fluid-structure interaction problems from pressure stresses to total stresses through the inclusion of viscous stress components, thus increasing the fidelity and versatility of TOFSI as a design tool.To overcome the discontinuities in the derivatives of the velocity across Lagrangian elemental boundaries, we utilized the superconvergent patch recovery technique by Zienkiewicz and Zhu to smoothen the velocity derivatives at the elemental boundaries.Sensitivity analysis equations of the implemented technique are derived and a sensitivity verification check is performed using the complex-step derivative approximation method.The key conclusions from this work are: • TOFSI problems using density-based methods are highly sensitive to the selection of the interpolation parameters in general and the Brinkman penalization interpolation parameter p α in specific.
• Depending on the flow characteristics, pressure vs. total stress force coupling may show a differentiation in the optimized designs, where each design works well for its intended purpose.Pressure spikes at sharp corners are a clear example.
• With efficient coding practices and in the absence of mesh deformations, the increase in computational costs when using total stress coupling as opposed to pressure coupling in TOFSI can be kept below 10%, which is a negligible factor when deciding the use of either coupling methods.With mesh deformations taken into account, the increase in computational costs is expected to be higher.
• The effect of the Reynolds number Re on its own cannot be isolated from its constituents.We have demonstrated with examples that controlling Re using the dynamic viscosity µ is entirely different from controlling Re using the characteristic velocity v max .We conjecture that the main reason behind this issue is that the forces experienced by the structure are not only related to Re but rather to its constituents (i.e.characteristic velocity, characteristic length, density, and dynamic viscosity).
A Numerical Verification of the Brinkman-Penalized Fluid-Structure Interaction Formulation In this appendix, we verify numerically the Brinkman-penalized, unified domain FSI formulation through a comparison with a body-fitted, segregated domain FSI formulation.To provide an additional layer of verification of our proprietary code, the body-fitted results are extracted from COMSOL Multiphysics [40].
The first case is an initial design of the column in a channel problem described in Sec. 4 assuming pressure coupling.The design domain is set to zero (i.e.ρ = 0 in Ω d ) while keeping the non-design column as solid (i.e.ρ = 1 in Ω nd ).To maintain the infinitesimal deformation assumption at this reduced material state, the elastic modulus is raised to 1e + 7 Pa.In COMSOL Multiphysics, we implement the same mesh and finite element type as our proprietary code, and the fluid-structure coupling is set to pressure forces only through modifying the weak expression under the fluid-structure interaction node.The results to be compared are the maximum and minimum x and y velocities, the maximum and minimum pressures, the total compliance of the structure, and the x and y displacements of a point 'A' located at the top left corner of the non-design column such that {x A , y A } = {0.975m, 0.5 m}.We also compare the fluid velocity field with streamlines and the pressure field with contour lines.As can be seen in Table 5, there is good agreement between the unified and segregated formulations for the parameters of interest.Similarly, the velocity and pressure fields show good agreement as shown in Fig. 13.
As for the total stress coupling, the forces experienced by the structure on the fluid-structure interface in the segregated domain formulation is from one side only, hence the discontinuity issue is not as severe as in the unified domain formulation implemented in this work.To the best of our knowledge, commercial software packages typically don't treat these discontinuities by default.In addition, discontinuity treatment using different techniques are not directly comparable.Hence, a similar verification with the total stress coupling case doesn't offer any value.
The second case is an optimized design of the column in a channel test problem using pressure coupling at Re = 1; namely Fig. 6c.To generate the geometry needed for COMSOL Multiphysics, the optimized design is traced manually in an external CAD software then imported into COMSOL Multiphysics.The fluidstructure coupling is set to pressure similar to the first case.An unavoidable byproduct of using the unified domain formulation, is the uncontrolled pressure inside the enclosed cavities.In the segregated domain formulation, this internal pressure doesn't typically occur.To generate the same boundary conditions in COMSOL Multiphysics, the fluid flow is first solved using the Brinkman-penalized Navier-Stokes to calculate these internal pressure values.This is performed by appending the fluid flow node with a volume force (i.e.−α max •v 1 in x and −α max •v 2 in y).Once these internal pressures are calculated in COMSOL, the volume force node is disabled and the entire problem is solved as a fluid-structure interaction problem in the segregated domain formulation where the calculated pressures are applied as boundary loads at the enclosed cavities.The calculated compliance using the segregated formulation is 0.01251 J against 0.01137 J for the unified domain formulation resulting in an error of 9.08%.There are mainly two reasons for this discrepancy; (i) the effect of the smooth interfaces implemented in COMSOL vs the original jagged interface in the unified domain solution, and (ii) the effect of some intermediate density elements on the outside interface that strongly affect the pressure field resolved.We believe both effects can be alleviated if some sort of adaptive mesh refinement is implemented during the optimization process, cf.[41,42].

Figure 2 :
Figure 2: Effect of the superconvergent patch recovery technique on velocity derivatives.

Figure 3 :
Figure 3: Location and number designations of the finite elements selected for sensitivity verification.

p 6 Figure 4 :
Figure 4: A comparison of two Brinkman penalization interpolation functions.Solid lines represent the interpolation function used by Yoon (Eq.33 in [18, p. 607]) while dashed lines represent Eq. 7 in this work.A unity α max and zero α min are used for demonstration purposes.

Figure 6 :
Figure 6: Effect of p E vs p Υ on optimized designs with p α = 2.5e − 6 at Re = 1 with pressure coupling.At iterations {21, 41, 61, 81}, p E is set to {1.5, 2.0, 3.0, 4.0} while p Υ is calculated for different δ E|Υ as shown in the captions.f is the objective function while DM is the discreteness measure.

Figure 7 :
Figure 7: Effect of p E vs p Υ on optimized designs with p α = 2.5e − 6 at Re = 1 with total stress coupling.At iterations {21, 41, 61, 81}, p E is set to {1.5, 2.0, 3.0, 4.0} while p Υ is calculated for different δ E|Υ as shown in the captions.f is the objective function while DM is the discreteness measure.

Figure 10 :
Figure 10: Effect of using v max to change the Reynolds number instead of µ.

Figure 12 :
Figure 12: Sensitivity of the column in a channel test problem at the second iteration for Re = 1 with p α = 2.5e − 6.

( a )
Velocity field with streamlines extracted from COMSOL Multiphysics.(b) Velocity field with streamlines extracted from our proprietary code.(c) Pressure field with contour lines extracted from COMSOL Multiphysics.(d) Pressure field with contour lines extracted from our proprietary code.

Figure 13 :
Figure 13: Numerical verification of the Brinkman-penalized FSI formulation for the initial design of the column in a channel problem.Velocity fields with streamlines and pressure fields and contour lines are shown for comparison.

Table 1 :
Verification of analytically-derived sensitivities using the complex-step derivative approximation method.
approximation method is implemented by introducing a small imaginary step ih to a particular design variable ρ i at a time, and calculating the approximate derivative of the objective function f w.r.t. that particular design variable ρ i as follows:

Table 2 :
Cross-check for the optimized designs using pressure coupling against Re.

Table 3 :
Cross-check for the optimized designs using total stress coupling against Re.

Table 4 :
Cross-check for the optimized designs against the type of force coupling.

Table 5 :
Numerical verification of the Brinkman-penalized FSI formulation for the initial design of the column in a channel problem.Some parameters of interest are selected for comparison.