Adjoint-based Shape Optimization for the Minimization of Flow-induced Hemolysis in Biomedical Applications
aa r X i v : . [ phy s i c s . f l u - dyn ] J a n Adjoint-based Shape Optimization for the Minimization of Flow-inducedHemolysis in Biomedical Applications
Georgios Bletsos , Niklas Kühl and Thomas RungHamburg University of Technology, Institute for Fluid Dynamics and Ship Theory, AmSchwarzenberg-Campus 4, D-21075 Hamburg, GermanyJanuary, 2021 Abstract
This paper reports on the derivation and implementation of a shape optimization pro-cedure for the minimization of hemolysis induction in biomedical devices. Hemolysis isa blood damaging phenomenon that may occur in mechanical blood-processing applica-tions where large velocity gradients are found. An increased level of damaged blood canlead to deterioration of the immune system and quality of life. It is, thus, important tominimize flow-induced hemolysis by improving the design of next-generation biomedicalmachinery. Emphasis is given to the formulation of a continuous adjoint complementto a power-law hemolysis prediction model dedicated to efficiently identifying the shapesensitivity to hemolysis. The computational approach is verified against the analyticalsolutions of a benchmark problem, and computed sensitivity derivatives are validated bya finite differences study on a generic 2D stenosed geometry. The application includedaddresses a 3D ducted geometry which features typical characteristics of biomedicaldevices. An optimized shape, leading to a potential improvement in hemolysis induc-tion up to 22%, is identified. It is shown, that the improvement persists for different,literature-reported hemolysis-evaluation parameters.
Keywords : Computational fluid dynamics (CFD); Adjoint-based shape optimization;Biomedical design; Hemolysis minimization [email protected] Introduction
The ever-growing advances in medicine, engineering and material science have led tothe development of biomedical devices, such as blood pumps, which allow long-termpatient care and significantly improve quality of life. Despite all the advances, a criticaltask for the design and development of such devices [33, 38] is still the minimization ofshear-induced blood damage (i.e. hemolysis) to guarantee good bio-compatibility.In the context of this manuscript, hemolysis refers to the mechanical damage of redblood cells due to excessively high stress induced by peculiarities of the blood flow. Itcan lead to hemoglobinemia, which plays a significant role in the pathogenesis of sepsis,and to increased risk of infection due to its inhibitory effects on the innate immunesystem [4]. Hemolysis induction is encountered in many biomedical devices, where largevelocity gradients are found [3], as well as in vivo conditions when vessels deliveringblood are kinked or stenosed [9]. The shape of the respective artificial devices or vesselsis believed to play a crucial part in the induction of blood damage due to its decisivefluid dynamic role. A computationally efficient shape optimization framework would betherefore appreciated for the development of next-generation biomedical machinery.The success of a numerical optimization process partially depends on the accuracy ofthe hemolysis-prediction model. While the study of shear-induced hemolysis has beenof interest for many experimental, in vitro conditions for quite some time [6, 12, 39], it isalso becoming increasingly important in silico conditions [30, 37]. Computational fluiddynamics (CFD) offer the possibility to predict blood damage in a purely numericalmanner and employ the damage prediction model as a cost function in an optimizationframework. A variety of such damage prediction models have been proposed [37]. Mostof them relate hemolysis with the magnitude of the shear stress and an exposure periodusing a variation of a power-law formulation, initially suggested by Giersiepen et al. [6].Despite the remarkable progress, a numerical model able to satisfactorily predict blooddamage in a variety of flows has not yet been established [9]. The present manuscriptdoes not aim to advocate the merits of a specific hemolysis evaluation or hemolysisprediction model, but is mainly concerned with the adaptation and integration of aclassical model into an adjoint shape optimization framework. Employing an Eulerianhemolysis-prediction model alongside a common CFD solver, we formulate a Partial-Differential-Equation (PDE) constrained optimization problem, which targets to min-imize flow-induced hemolysis by improving the shape. In the context of this paper,hemolysis will be referred to as objective while the shape as control . In CFD-basedoptimization, multiple ways, ranging from stochastic [1, 25] to deterministic [2, 34] opti-mization methods, could be followed. The present research is concerned with the efficientcomputation of the derivative of the objective with respect to (w.r.t) the control. Theaforementioned derivative is subsequently used by a deterministic gradient-based steep-est descent method, which drives the controlled shape towards an improved state. Tothat extent, the continuous adjoint method is studied. The adjoint method has been,increasingly, receiving attention in terms of CFD-based optimization [11, 14, 22, 23, 31]2ince the pioneering works of Pironneau [26] and Jameson [15], due to its superior com-putational efficiency. In specific, the attractiveness of the method lies on the fact thatthe computation of the derivative is independent from the size of the control. For fur-ther details on the merits and drawbacks of the adjoint approach, the interested readeris referred to [7, 27].The remainder of this paper is organised as follows: Section 2 presents the mathemat-ical model of the primal and adjoint problem. The same section also reports on theemployed numerical procedure and utilized boundary conditions (BCs). In Section 3, abenchmark case is studied to verify the code reliability. Moreover, the accuracy of thecomputed sensitivity derivative is assessed in the context of a Finite Differences (FD)study. Subsequently, in Section 4, the model is applied on a 3D geometry. Results for afull shape optimization process are presented and their dependency on parameters of thenonlinear hemolysis evaluation model is discussed. The paper closes with conclusionsand outlines further research in Section 5.Within this publication, Einstein’s summation convention is used for repeated lower-caseLatin subscripts. Vectors and tensors are defined with reference to Cartesian coordi-nates.
This section is dedicated to the formulation of the primal (physical) and adjoint (dual)problem. The coupling of the hemolysis-prediction model with the primal flow equationsis discussed. A presentation of the newly developed adjoint model then follows withdetailed discussions on adjoint-hemolysis specific points of interest.
Throughout this paper, blood is treated as an incompressible, Newtonian fluid. Theassumption of non-Newtonian behaviour has been shown to satisfactorily predict theflow of blood in tubes with diameter ranging from 130 to 1000 µm [24]. However, allthe applications considered herein refer to ducted systems with a larger diameter andthus the Newtonian assumption is preferred. Furthermore, we assume all applications inlaminar and steady-state conditions. The blood flow in a domain (Ω) follows from theNavier-Stokes (NS) equations for the conservation of volume and momentum, viz. R p = − ∂u j ∂x j = 0 , (1) R ui = ρu j ∂u i ∂x j − ∂∂x j [2 µS ij − pδ ij ] = 0 . (2)3ere u i , p , S ij = 0 . ∂u i /∂x j + ∂u j /∂x i ), δ ij , ρ and µ refer to the components of thefluid velocity vector, static pressure, components of the strain-rate tensor, Kroneckerdelta components as well as the fluid density and dynamic viscosity, respectively.In the framework of CFD-based hemolysis analysis, a typical approach is to utilize aprediction or evaluation model in a one-way coupling with the CFD solver [37]. Thisimplies that the hemolysis model receives information from the NS equations (1, 2) withno retro-action on the employed fluid properties ( ρ, µ ). The hemolysis prediction modelconsidered in this study originates from the power-law equation, first introduced byGiersiepen et al. [6], and reads H ( τ , t ) = Cτ α t β . (3)The hemolysis index H denotes a measure for the released hemoglobin to the totalhemoglobin within the red blood cell. It is governed by a scalar stress representation τ and an exposure time t , which represents the duration on which the red blood cell isexposed to the stress. The remaining constants ( C , α , β ) are introduced to fit experi-mental data. Table 1 summarizes the parameter sets utilized in this study. The secondNotation Values ( C, α, β ) Type of Blood Max. Stress (Pa)GW (3 . · − , 2 . . . · − , 1 . . . · − , 1 . . τ ij = 2 µS ij is frequently used to determine the scalar stressrepresentation τ in combination with a user specific parameter k = 1 , , τ = p − kI τ with I τ = 12 (cid:2)(cid:0) tr ( τ ij ) − tr ( τ ij ) (cid:1)(cid:3) . (4)Due to tr ( τ ij ) = tr (2 µS ij ) ∼ ∂u i /∂x i , the first contribution to I τ in (4) vanishes forincompressible Newtonian fluids, which yields τ = (cid:0) kµ S ij S ji (cid:1) := τ ∗ . (5)To incorporate Eqn. (3) into an adjoint optimization framework, a PDE-based hemolysisprediction is advantageous. To this end, a linearization strategy is used [5] for (3), i.e. H β = (cid:16) C β τ αβ (cid:17) t. (6)Substituting H L = H β , the sought material derivative expression follows from (6) ρ ∂H L ∂t + ρ u j ∂H L ∂x j = ρ C β τ αβ . (7)4ssuming steady-state, a modification of the residual form of (7) reads R H = ρ u j ∂H L ∂x j − ρ C β τ αβ (1 − H L ) = 0 in Ω . (8)Note that (8) is enhanced by (1 − H L ) to ensure H L <
1. Equation (8) will be referredto as (primal) hemolysis equation throughout this manuscript and serves to formulatethe objective functional of the optimization. The merits of expression (8) refer to itsstraightforward derivation from the original power-law formulation (3) and its suitabilityfor showcasing the adjoint method in the context of blood damage. While it is true thathemolysis does not occur for small values of τ < τ th , where τ th corresponds to somethreshold value (which might also be related to the exposure time), this aspect is delib-erately ignored in the present effort. Furthermore, the footing of the hemolysis modelassumes a homogeneous stress representation. Even though this is a strong assumptionfor most practical applications, the predictions have been shown to be sufficiently reli-able when relatively comparing different geometries [32]. In the present context of shapeoptimization, it is thus fair to assume that the qualitative information on a reduction ofthe induced hemolysis is reliable.Equation (8) serves for the formulation of an adjoint complement in Sec. 2.2.1. Fur-thermore, the optimization requires a scalar objective functional. An integral indexsignifying the level of flow-induced hemolysis can be computed from the domain integral R ∇ j ( ρu j H L ) d Ω, or a corresponding boundary integral. Assuming zero inflow fluxes of H L , no-slip velocity on the walls and solenoidal flow fields, only the outlet fluxes remain,which are usually normalized by the respective mass flux H index = R Γ out Hρu j n j d Γ R Γ out ρu j n j d Γ = R Γ out H βL ρu j n j d Γ R Γ out ρu j n j d Γ . (9)Note that the denominator of Eqn. (9) is constant in steady flows of incompressible fluids.Hence, the numerator of (9) is considered as an appropriate objective or cost-functionalin this manuscript. This subsection outlines the formulation of the adjoint problem at hand which is initiallyexpressed as an optimization problem constrained by a set of PDE, namely R ( y, c ) = 0,viz. min c ∈ C ad J ( y, c ) s.t. R ( y, c ) = 0 . (10)Here J is the objective or cost functional, y is the vector of primal state variables, whichconsists of p , u i and H L , and c is the control parameter which we need to find in the setof admissible states C ad . In this study, c refers to the shape of the structure in whichblood flows. It is thus a subset of ∂ Ω or Γ -the latter notation is preferred throughout5his manuscript- and can be further restricted by considering only specific sections ofthe structure, namely Γ D ⊂ Γ, where the index D denotes design.Based on (9) the objective functional under investigation can be written as J = Z Γ out H βL ρu i n i d Γ = Z Γ out j Γ d Γ , (11)while the set of PDE, R ( y, c ) = 0, consists of Eqns. (1), (2) and (8). Having formulatedthe optimization problem as a constrained problem in Eqn. (10), the Lagrange principleto eliminate the constraints by employing appropriate Lagrange multipliers is utilized[29]. In this context, a continuous Lagrange functional is defined as L := J + Z Ω h ˆ pR p + ˆ u i R ui + ˆ hR H i d Ω , (12)where ˆ p , ˆ u i and ˆ h are Lagrange multipliers that are referred to as adjoint pressure,components of the adjoint velocity vector and adjoint hemolysis, respectively. Theirdimensions [ˆ u i ] = [ J ] / ([ R u i ] m ), [ˆ p ] = [ J ] / ([ R p ] m ) and [ˆ h ] = [ J ] / ([ R H ] m ) depend onthe underlying cost functional, where [ J ] = ([ j Γ ] m ) represents the units of the integralboundary-based objective. Since, from a physical point of view, equations R p , R ui and R H are required to be zero, it is apparent that L is equal to J . The optimization problem(10) is, thus, equivalent to min L (ˆ y, y, c ), where ˆ y is the vector of adjoint variables and y is no longer constrained. For a minimum of (12) the total variation δ L vanishes δ L = δ ˆ y L · δ ˆ y | {z } FD1 + δ y L · δy | {z } FD2 + δ c L · δc | {z } FD3 = 0 . (13)The terms FD1, FD2 and FD3 denote the variation of the Lagrangian in the directionof adjoint ( δ ˆ y ) and primal ( δy ) state as well as the control ( δc ), respectively. They canbe computed by utilising the functional derivative into the respective direction. FD1leads to the known set of primal equations that need to be satisfied for every controlstate. FD2 yields the accompanying adjoint equations. Finally, FD3 gives rise to asensitivity derivative that offers information of the objective w.r.t. the control. Eqn.(13) is only satisfied when a global or local minimum is reached. Hence, the followingthree optimality conditions are obtained for vanishing FD2 δ p L · δp = 0 ∀ δp , δ u i L · δu i = 0 ∀ δu i , δ H L L · δH L = 0 ∀ δH L . (14)A detailed overview on PDE constraint optimization can be found in [13]. Each optimality condition (14) leads to a PDE that needs to be satisfied in order forFD2 to vanish. The essence of the adjoint method is to identify an adjoint state ˆ y that satisfies (14), so that the total variation δ L depends only on the variation of the6ontrol ( δ c L · δc ). For the purpose of this paper it is deemed beneficial to split (12) intosub-integrals consisting of isolated terms of the primal equations, viz. Z Ω ˆ p (cid:16) − ∂u j ∂x j (cid:17) d Ω | {z } I + (cid:20) Z Ω ˆ u i (cid:16) ρu j ∂u i ∂x j (cid:17) d Ω | {z } I + Z Ω ˆ u i (cid:16) − ∂∂x j (2 µS ij ) (cid:17) d Ω | {z } I + Z Ω ˆ u i (cid:16) ∂p∂x i (cid:17) d Ω | {z } I (cid:21) + (cid:20) Z Ω ˆ h (cid:16) u j ∂H L ∂x j (cid:17) d Ω | {z } I + Z Ω ˆ h (cid:16) − C β τ αβ (1 − H L ) (cid:17) d Ω | {z } I (cid:21) . The computation of the functional derivatives (14) requires the computation of the deriva-tives of each sub-integral I k , with k = 1 , ...,
6, as well as the derivatives of the objectivefunctional (11). The latter involves δ p J · δp = 0 and δ u i J · δu i = Z Γ out δu i (cid:16) ρH βL n i (cid:17) d Γ , δ H L J · δH L = Z Γ out δH L (cid:16) ρβH ( β − L u i n i (cid:17) d Γ . (15)Note that the objective functional and its directional derivatives only exist at the outletof the domain, which is not a design surface (Γ D ). The contribution of the objectivefunctional to the adjoint equations thus appears only in the outlet boundary conditions.In an analogy to primal duct flows which are driven by a difference between the inletand the outlet pressure, the dual flow is driven by the relative difference of the hemolysisindex.The calculation of the variations δ y I k · δy , with k = 1 , .., I and I represents a core contribution of this paper. For δ y I · δy , we can utilize Gauss’sdivergence theorem and obtain δ p I · δp = 0 (16) δ u i I · δu i = Z Ω δu i (cid:16) ˆ h ∂H L ∂x i (cid:17) d Ω = Z Γ ( δu i n i ) H L ˆ h d Γ − Z Ω δu i (cid:16) H L ∂ ˆ h∂x i (cid:17) d Ω (17) δ H L I · δH L = Z Ω ˆ hu j ∂δH L ∂x j d Ω = Z Γ δH L ( u j n j )ˆ h d Γ − Z Ω δH L (cid:16) u j ∂ ˆ h∂x j (cid:17) d Ω . (18)Although Eqn. (17) doesn’t necessarily require integration by parts to be included inthe adjoint equations, the respective volume integral is expanded in the context of thiswork for computational purposes.Following the same strategy for I , yields δ p I · δp = 0. As for δ u i I · δu i , we first define M = C β (1 − H L ) as a scalar quantity not subjected to any variation w.r.t the velocity.Using the scalar stress expression (5), the integral I is expressed as I = − R Ω ˆ hM τ ∗ α β d Ω7nd its linearized variation in the velocity direction reads δ u i I · δu i = − Z Ω ˆ hM α β τ ∗ ( α β − ( δ u i τ ∗ · δu i ) d Ω= − Z Ω ˆ h (cid:16) M α β τ ∗ ( α β − kµ S ij (cid:17) δS ij d Ω , (19)where δS ij = 0 . ∂δu i /∂x j + ∂δu j /∂x i ). By defining B ij = 2 M kµ αβ τ ∗ ( α β − S ij , acompact form of (19) is further expanded via integration by parts as δ u i I · δu i = − Z Ω ˆ h B ij (cid:16) ∂δu i ∂x j + ∂δu j ∂x i (cid:17) d Ω = − Z Ω ˆ hB ij (cid:16) ∂δu i ∂x j (cid:17) d Ω= − Z Γ δu i (cid:16) ˆ hB ij n j (cid:17) d Γ + Z Ω δu i ∂ (ˆ hB ij ) ∂x j ! d Ω . (20)Finally, δ H L I · δH L is calculated as δ H L I · δH L = Z Ω δH L (cid:16) ˆ hC β τ αβ (cid:17) d Ω . (21)Having expressed the variation of J and I k in terms of several boundary and volumeintegrals, we can superpose the variation of the Lagrange functional as the sum of oneboundary and volume integral, viz. δ y L · δy = Z Γ δy ( · ) d Γ + Z Ω δy ( · ) d Ω . (22)The adjoint equations as well as their corresponding boundary conditions arise by de-manding that the integrands in Eqn. (22) vanish for any test function δy . Therefore,following the aforementioned methodology for the optimality conditions (14), the result-ing constraints in the interior of the domain (Ω) readˆ R p = − ∂ ˆ u i ∂x i = 0 , (23)ˆ R ui = ρ ˆ u j ∂u j ∂x i − ρu j ∂ ˆ u i ∂x j − ∂∂x j (cid:16) µ ˆ S ij − ˆ pδ ij (cid:17) + ∂ (ˆ hB ij ) ∂x j − H L ∂ ˆ h∂x i = 0 , (24)ˆ R H = u j ∂ ˆ h∂x j − C β τ αβ ˆ h = 0 . (25)Equations (23) as well as (24) represent the adjoint companions to the continuity andmomentum equation, respectively. In this sense, Eqn. (25) is referred to as adjointhemolysis equation throughout this manuscript. The adjoint momentum equation isenhanced with the last two terms on the left-hand-side (LHS) of the equation thatinclude contributions from adjoint and primal hemolysis. This is due to the fact thaton the primal side of the system, the velocity appears twice in the hemolysis equation,8nce in the advection and once in the form of a gradient in the source term. Both termsact similar to a pressure term in (24) and thereby drive the adjoint flow field.Furthermore, it is also worth saying that in contrast to its primal counterpart, theadjoint hemolysis equation doesn’t require any information from the solution of theadjoint continuity and momentum equations. Once again, we have a one-way couplingbut this time the direction of the coupling and thereby the algorithmic order of sequenceis reversed. Practically this means that the adjoint hemolysis equation could be solvedafter the solution of the primal equations and prior to the solution of the other twoadjoint Eqns. (23) & (24).The adjoint system (23-25) is an extension of a classical adjoint (incompressible) NSsystem, which corresponds to (23-24) excluding the two hemolysis terms. This is similarto the adjoint complement of the Volume of Fluid (VoF) method [16, 18], used for multi-phase applications, and enables a straightforward implementation, provided that anadjoint multi-phase solver exists. Although the adjoint equations are per definition linear,several cross coupling terms might introduce a severe stiffness to the densely coupledPDE system. A viable workaround to stabilize the iterative procedure refers to theintroduction of additional diffusive terms or adjoint pseudo time-stepping, cf. [8,16,18].The adjoint BCs arise by fulfilling the requirement of eliminating the surface integral of(22). By expanding this necessity for every primal state variable the BCs readˆ u i n i = 0 along Γ , − ˆ pn i + ρ ˆ u i u j n j + µ ∂ ˆ u i ∂x j n j + H L ˆ hn i − ˆ hB ij n j = 0 along Γ \ Γ out , − ˆ pn i + ρ ˆ u i u j n j + µ ∂ ˆ u i ∂x j n j + H L ˆ hn i − ˆ hB ij n j + ρH βL n i = 0 along Γ out , − ˆ h ( u j n j ) = 0 along Γ \ Γ out , ˆ h + ρβH ( β − L = 0 along Γ out . (26) Ensuring a vanishing residual of the primal and adjoint system of equations leads tovanishing FD1 and FD2 terms of (13). The FD3 term allows for the computation ofa surface sensitivity which is used by a gradient-based algorithm to drive the shapetowards an improved state. Following (12), FD3 is expressed as δ c L · δc = δ c J · δc + Z Ω ˆ y (cid:16) δ c R · δc (cid:17) d Ω . (27)Since the primal flow equations (1, 2, 8) are also fulfilled under the total variation, oneobtains δR ( y, c ) = δ c R · δc + δ y R · δy = 0 . (28)9herefore, the variation w.r.t. the control transforms to δ c R · δc = − δ y R · δy. (29)Furthermore, based on the formulation of our objective functional J , which is definedonly on the outlet, δ c J · δc = 0 since c ≡ Γ D and Γ D ∩ Γ out = ∅ . The sensitivity derivativecan thus be written as δ c L · δc = − Z Ω ˆ y (cid:16) δ y R · δy (cid:17) d Ω . (30)The right-hand-side (RHS) of (30) is developed further and by taking into considerationthe BCs of the problem as well as a Taylor expansion of the velocity w.r.t a perturbationof the design wall (see App. B or [20]), the sensitivity to be computed reads S L = Z Γ D − (cid:20) µ ∂u i ( t ) ∂n ∂ ˆ u i ( t ) ∂n (cid:21) d Γ , (31)where u i ( t ) , ˆ u i ( t ) denote the part of the primal and adjoint velocity tangential to thesurface, respectively and n denotes the normal to the surface. Interestingly, the shapederivative is not directly affected by the adjoint or primal hemolysis fields. Nevertheless,the primal and adjoint hemolysis parameters drive the adjoint flow field, cf. (24), andthereby propagate into the shape derivative through the adjoint velocity gradient. In what concerns the application studied in this work, the flow boundary Γ consistsof three parts, namely the inlet (Γ in ), outlet (Γ out ) and wall (Γ W ). We split the wallboundary into two portions, Γ D and Γ B = Γ W \ Γ D , where the first one involves the partsunder design while the latter one is bounded to the initial configuration. The completeboundary domain is thus described as Γ = Γ in ∪ Γ out ∪ Γ D ∪ Γ B .While the BCs of the primal system of equations are selected based on the physicalproperties of the problem, the adjoint BCs are deduced based on the necessity of van-ishing the surface integrals of Eqn. (22). In many boundary patches this necessity isinherently satisfied due to the primal BCs, though in general it reduces to satisfyingEqns. (26). Furthermore, it is worth mentioning that the BCs strongly depend on theobjective functional, if that is defined in a subset of Γ. The complete set of BCs for theprimal and adjoint problem are summarized in Tab. 2.10oundary Patch u i p H L ˆ u i ˆ p ˆ h Γ in u i = u ini ∂p∂n = 0 H L = 0 ˆ u i = 0 ∂ ˆ p∂n = 0 ∂ ˆ h∂n = 0Γ out ∂u i ∂n = 0 ∂p∂n = 0 ∂H L ∂n = 0 ˆ u i = 0 2nd Eqn. (26) 5th Eqn. (26)Γ W u i = 0 ∂p∂n = 0 ∂H L ∂n = 0 ˆ u i = 0 ∂ ˆ p∂n = 0 ∂ ˆ h∂n = 0Table 2: Boundary conditions for closing the primal and adjoint equation system. The numerical procedure for the solution of the primal and adjoint system is based uponthe Finite Volume Method (FVM) employed by FreSCo + [28]. Analogue to the use ofintegration-by-parts in deriving the continuous adjoint equations [19,21], summation-by-parts is employed to derive the building blocks of the discrete adjoint expressions. Adetailed derivation of this hybrid adjoint approach can be found in [16, 18, 31]. The lasttwo terms of the adjoint momentum equation, involving hemolysis contributions, areadded explicitly to the RHS. The segregated algorithm uses a cell-centered, collocatedstorage arrangement for all transport properties. The implicit numerical approximationis second order accurate in space and supports polyhedral cells. Both, the primal andadjoint pressure-velocity coupling is based on a SIMPLE method and possible paral-lelization is realized by means of a domain decomposition approach [35, 36]. The complete shape-optimization framework is summarized in Fig. 1. The efficiency ofthe method is illustrated by taking into consideration the required computational budgetfor one complete optimization process based on the presented flowchart. Assuming thatthe cost for solving the primal problem is equal to 1 equivalent flow solution (EFS) thenthe cost for solving the adjoint one is approximately also equal to 1 EFS. The final twoprocesses (represented by rectangles on the flowchart) have a practically negligible costcompared to the EFS unit. The complete optimization cost is thus equal to i max times2 EFS, regardless of the size of the control which for our case is proportional to thenumber of discretized nodes on Γ D . In the context of this study, the steepest descentmethod is considered to advance the shape.The computed surface sensitivity S L might possibly be rough. We thus employ theLaplace-Beltrami [17] metric to extract a smooth gradient G L i through a numericalapproximation of G L i − λ ∆ Γ G L i = S L n i , (32)where ∆ Γ refers to the Laplace-Beltrami operator (∆ Γ = ∆ − ∆ n ), λ corresponds to auser-defined control of smoothing and n i denotes the normal vector at each face of Γ.11tartInitial shape c ,Solution settings (fluidproperties, numer-ical schemes, e.t.c)Solve primal problem in c i ,Eqns. (1), (2),(8) −→ Acquire y .Solve adjoint problemin c i based on known y ,Eqns. (23), (24),(25) −→ Acquire ˆ y .Compute surface sensitivitybased on known y and ˆ y ,Eqn. (31) −→ Acquire S L .Update shape to c i +1 usingthe steepest descent method. Is i < i max ? EndYES NOFigure 1: Flowchart of adjoint shape optimization framework for one design candidate i .Once the smooth gradient field is available, the displacement field d i is computed fromEqn. (33), ∂∂x j (cid:20) q ∂d i ∂x j (cid:21) = 0 in Ω , d i = G L i along Γ D , d i = 0 along Γ \ Γ D . (33)where q = 1 / ( W D + ǫ ), W D is the wall normal distance and ǫ = 10 − .12 Verification and Validation Studies
This section verifies the implementation of the primal and adjoint hemolysis systemfor a benchmark problem, as suggested by the Food & Drug Administration (FDA) [10].Analytical solutions are derived and compared with numerical predictions. Subsequently,a finite difference study is conducted on a 2D geometry to validate the sensitivity.
The benchmark problem refers to a fully developed pipe flow, considered on a 3D mesh.For brevity reasons, the derivation of the analytical primal hemolysis solution is skippedand the reader is referred to [10]. The solution of the primal flow reads u z ( r ) = U max (cid:18) − (cid:16) rR (cid:17) (cid:19) , H ( r, z ) = C τ α z β u βz , τ ( r ) = 2 µU max rR , (34)where r corresponds to the radial direction, R marks the pipe radius, U max refers tothe centerline velocity, z to the axial coordinate and the direction of the primary flowvelocity u z .The considered verification case is sketched in Fig. 2. The entrance and exit planes ofthe pipe refer to z = 0 and z max = 2 m . The pipe radius reads R = 0 . m . The fullydeveloped axial velocity profile of the laminar flow utilizes U max = 10 m/s . The fluidproperties are selected based on a prescribed Reynolds number of Re = 2 U R/ν = 2000,where U and ν refer to the bulk velocity and kinematic viscosity, respectively. Thethree-dimensional geometry is discretized with approximately 75k control volumes ona structured grid. A cross-section of the computational grid employed is shown in Fig.3.Due to the decoupled nature of the adjoint hemolysis equation (25), an analytical solutioncan be stated using the analytical primal flow solution (34), viz. u z d ˆ hdz − C β τ αβ ˆ h = 0 . (35)If we abbreviate Λ = C β τ αβ , then the solution of (35) readsˆ h = K e Λ uz z , (36)The integration constant K is calculated based on the BCs of the adjoint hemolysisequation which results from satisfying the final expression in (26) and demands ˆ h | out = − ρβH β − L | out . Finally, the adjoint hemolysis field in a fully developed pipe flow readsˆ h = − ρβH β − L (cid:12)(cid:12)(cid:12) out e Λ uz ( z − z max ) , (37)13 = 0.5 m z max = 2 mU max = 10 m/s Figure 2: Sketch of the pipe flow bench-mark. Figure 3: Cross-section of the pipe meshcoloured with the axial velocity.where H L | out is calculated by setting z = z max in (34) together with H βL = H .Figure 4 compares the computed primal hemolysis solution against analytical valuesas calculated by Eqn. (34). The comparison is realized for all three sets of hemolysismodel parameters outlined in Tab. 1. As can be seen, all computed values are fittingthe analytical solutions to a satisfying degree. Figure 5 compares the computed adjoint -14 -13 -12 -11 -10 -9 -8
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 H e m o l ys i s ( H ) r (m) Analytical, GWSimulation, GWAnalytical, HOSimulation, HOAnalytical, ZTSimulation, ZT Figure 4: Primal hemolysis profiles at z =1 . m for three different parame-ter sets ( C , α , β ). The notationfollows from Tab. 1. Continuouslines correspond to computationswhile symbols to analytical values.
100 1000 10000 100000 1x10
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 A b s o l u t e A d j o i n t H e m o l ys i s ( | ˆ h | ) r (m)Analytical, GW, z = 0.5mAnalytical, GW, z = 1.5mSimulation, GWAnalytical, HO, z = 0.5mAnalytical, HO, z = 1.5mSimulation, HOAnalytical, ZT, z = 0.5mAnalytical, ZT, z = 1.5mSimulation, ZT Figure 5: Magnitude of adjoint hemolysisprofiles for three different param-eter sets ( C , α , β ). Notation asin Fig. 4. Analytical results arereported for two longitudinal posi-tions, z = 0 . z = 1 . m .hemolysis field against the analytical solutions (37). Due to the nature of the BCs of theadjoint hemolysis equation, the field values are negative and thus the absolute values14re presented to support a log-scale of the ordinate. Predictions are nearly identicalto analytical solutions for all three set of parameters. Both the primal and adjointhemolysis solvers are thus considered verified.We would like to remark on the nature of the adjoint hemolysis equation. As can beseen in Fig. 5, the adjoint hemolysis profile changes only slightly for different axialpositions. The BC at the outlet, which reads ˆ h | out = − ρβH β − L | out , dominates thecomplete field since all three cases correspond to β < H L <<
1. To avoidnumerical errors that would arise for H L | out = 0, the BCs on the outlet is reformulatedto ˆ h | out = − ρβ ( H L + ǫ ) ( β − | out , where ǫ = 10 − . Based on the previous comments, onecould assume that applying the bulk outlet adjoint hemolysis profile to the whole fieldsuffices. However, due to the existence of the final two terms in the adjoint momentum,(cf. Eqn. (24)), this assumption would fall short on capturing the conceptual descriptionof the complete model. The goal of the primal-adjoint simulation is the computation of the surface sensitivity(shape derivative) through Eqn. (31). It is thus important to investigate the accuracyof the computed sensitivity. In order to realize this, the computed shape derivative iscompared against locally evaluated second order accurate finite differences [14,17], viz. δ c i J = [ J ( c i + ǫn i ) − J ( c i − ǫn i )]2 ǫ . (38)Here c i represent discrete points of the control, ǫ is the magnitude of the perturbationand n i is the normal vector at c i . In practise, the study is realized by deforming theboundary faces of the discretized geometry into their normal direction with a magnitudeequal to ǫ . The local boundary perturbations are then transported into the domainbased on (33). Figure 6 presents a schematic of the considered symmetric geometry. Inspecific, the design section (Γ D ) follows from y = A (cid:16) sin (cid:16) πL x − π (cid:17)(cid:17) , where A is theheight of the bump, L is the length of the design surface and L/A = 20. The domain isdiscretized with approximately 70k control volumes on a structured grid which is refinedin normal as well as tangential direction towards the wall and Γ D . A detail of the utilizedgrid along a part of the design section is shown in Fig. 7.The study is conducted for a parabolic inlet velocity profile with U max = 0 . m/s . Thefluid properties ( ρ, µ ) are assigned to unity to simplify the study. The results of the FDstudy are shown in Fig. 8 for two perturbation magnitudes, ǫ , at 20 selected discretepositions. The overall agreement between the computed adjoint shape derivative andthe FD results is very good. Furthermore, as shown in Fig. 9 the computed objectivefunctional, on the perturbed shapes, exhibits a linear behaviour w.r.t the perturbationsize.Overall, the study validates, on a preliminary basis, the accuracy of the adjoint method,as presented within this manuscript. 15 . Inlet 0.5 0.50.1
L AA = 0.005 Γ B Γ B Γ D Figure 6: Geometry for the 2D FD study.Blue lines correspond to the de-sign walls (Γ D ) while black linesindicate bounded walls (Γ B ). Alldimensions in (m). (Figure not inscale). Figure 7: Detail of the utilized grid for theFD study at the section x ∈ [0 . − . m , where x = 0 refersto the inlet and x = 1 . -25-20-15-10-5 0 5 0.5 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.6 S en s i t i v i t y D e r i v a t i v e ( [ J ] • m - ) L (m) Continuous adjointFD : ε /A = 10 -6 FD : ε /A = 10 -5 Figure 8: Comparison of sensitivity pre-dicted by Eqn. (31) (continuousline) and computed from FD, Eqn.(38), for ǫ/A = 10 − (circles) and ǫ/A = 10 − (triangles). -3-2-1 0 1 2 3 4-10 -1 1 10 ( J - J * ) • / J ( ε / A) • 10 ε / A = 10 -5 ε / A = 10 -6 Figure 9: Influence of perturbation size ǫ onFD results at x = 0 . m . J and J ∗ represent the objective func-tional of the unperturbed and per-turbed shapes, respectively. Having assessed the code implementation and the reliability of the computed sensitivity,the approach is put to the test by considering a complete optimization process on a3D geometry. The latter corresponds to a benchmark model, designed by a technicalcommittee to include flow phenomena related to blood damage in medical devices [30].16 .1 Investigated Configuration
The geometry under consideration shares characteristics of many blood-carrying medicaldevices. It includes sections where the flow is accelerating or decelerating, where re-circulation and significant variations in shear stress occur. All of these phenomena arebelieved to be related to blood damage in medical devices [30]. A merit of the adjointmethod is that it stems from the primal (physical) problem and does not operate asa black box. It is therefore capable of shedding light into areas of potentially elevatedblood-damage capacity. Areas on which an accumulated shape sensitivity is identifiedare expected to be physically relevant to the problem of hemolysis in general.Geometrical details of the model are presented in Fig. 10. The wall boundary is againsplit into 2 parts, Γ D and Γ B , as was done in the previous section. The inlet and outlettubes are considered to be bounded, since their shape is believed to be trivial to theblood damage capacity. The remaining structure is classified as Γ D and presented inblue. Inlet 80 22 40 80 Γ B Γ B Γ D d = 4 d = 12 Figure 10:
Top : Sketch of the investigated geometry, adapted from [30]. Presented inblue are the sections free for design. The remaining structure is bounded toits initial configuration. All dimensions in (mm).
Bottom : Detail of the computational grid on a longitudinal section of thegeometry near the design surface.The geometry is discretized with approximately 1 million control volumes on a butterfly-like structured grid. The mesh is gradually condensed near the walls, cf. Figs. 11 & 12,to adequately resolve relevant flow phenomena and also near the Γ D section (Fig. 10),to ensure an accurate computation of the sensitivity. As can be seen in Fig. 12, thegrid is additionally refined in the throat region to sufficiently capture the free shear flow,occurring by the jet exiting the throat.The fluid’s density and dynamic viscosity are set to 1056 kg/m and 3 . cP , respec-17ively, representing blood under physiological conditions. A fully developed laminaraxial velocity profile is prescribed at the inlet of the geometry, based on (34), so that theReynolds number at the throat reads Re T = 500. As regards the hemolysis-predictionmodel, the utilized set of parameters corresponds to GW (cf. Tab. 1) and we employ k = 1 for the computation of the scalar stress representation (4). The diffusion termsin the adjoint and primal momentum equations are discretized using the second orderaccurate central difference scheme while the convective terms are discretized through thehigher order Quadratic Upstream [Downstream] Interpolation of Convective Kinematics(QU(D)ICK) scheme.Figure 11: Cross-sectional view of the gridnear the inlet. Figure 12: Cross-sectional view of the gridnear the outlet. The optimization study was performed for 50 design candidates. The employed smooth-ing parameter value (32) reads λ = d T , where d T is the throat diameter. To avoid largedeformations of the geometry and the grid, the local displacement field d i computedfrom Eqn. (33) was scaled to a user-defined constant maximum value ˜ d max , viz.˜ d i = d i max ( d i ) ˜ d max . (39)For the study at hand ˜ d max /d T = 0 . · − .Figure 13 shows the displacement magnitude on the design surfaces of the structure athand, after the solution of the first primal and adjoint problem. As can be seen, thedisplacement is accumulated, almost entirely, at the sudden expansion of the geometry,18here the highest values of hemolysis occur. Furthermore, re-circulation is occurringwith relatively significant values of upstream mass fluxes after the expansion. This isbelieved to be one of the causes of hemolysis. Therefore, even though there is no directcontribution from the primal or adjont hemolysis to the shape sensitivity, the necessaryinformation from both fields are preserved through the adjoint velocity, which is directlyinfluenced by both H L and ˆ h .Figure 13: Displacement field as computed by Eqn. (39) after the first primal and adjointsimulations.The optimization history is shown in Fig. 15. It can be seen that the optimization startsconverging after approximately 40 shapes. The final shape results in a relative reductionof the objective functional by 22%. A comparison of the initial and the optimized shape isdisplayed in Fig. 14. The optimization algorithm, proceeded into widening and relativelysmoothing the sudden expansion of the geometry. This results in a smaller maximumvalue of the axial velocity as can be seen in Fig. 16. At the same time, nonetheless, theflow inside the bounded portions of the structure ( z ∈ [0 : 0 . m ) remains relativelyunchanged as shown on the same figure. This is important in terms of biomedicalapplications due to their ultimate goal of realizing a sensitive task. However, due tothe smoothing of the expansion zone, the velocity profiles are changed in the perturbedsection of the geometry (cf. Figs. 17 & 18). The velocity profile of the optimizedshape is smoothed near the wall region resulting in substantially lower shear stresses.Subsequently, due to the direct relation between shear stresses and hemolysis induction,the maximum values of hemolysis in the flow is reduced in the perturbed section of thestructure, which is also the area on which the maximum values are identified in theinitial shape. A direct comparison of the hemolysis profiles at two cross sections of the19 Symmetry line Y ( m ) Z (m)
Initial shapeOptimized shape
Figure 14: Projection of Γ D on a 2D plane.Red line illustrates the outline ofthe final optimized shape ( c )while black line corresponds tothe initial shape. The bottompart is symmetric w.r.t the sym-metry line. Axes not in scale. -25-20-15-10-5 0 0 5 10 15 20 25 30 35 40 45 50 ( J i - J ) / J x % Number of Shapes
Figure 15: Optimization convergence. J i de-notes the objective functional ascomputed for each shape c i and J denotes the same value for theinitial shape, c . Points illustratethe computed quantity at eachdiscrete shape. U z ( m / s ) Z (m) Initial shapeOptimized shape
Figure 16: Axial velocity along model’s centerline ( r = 0) for initial (black line) andoptimized (red line) shape. Coordinate Z = 0 corresponds to the model’sinlet.initial and final geometry is presented in Figs. 19 & 20.As described in the previous section, the optimization utilized a specific set of parametersfor the primal hemolysis model. It is interesting, therefore, to examine the performance ofthe optimized shape for other sets of hemolysis specific parameters. As shown in Fig. 21,for any choice of parameters mentioned in Tab 1, the optimal shape always outperformsthe initial one in terms of flow-induced hemolysis. In specific, the most significantimprovement occurs for the employed set of parameters during the optimization run,namely GW , while the value of k is relatively trivial w.r.t the improvement. This showsa robustness of the method in terms of user-defined values.20 U z ( m / s ) r (m) Initial shapeOptimized shape Figure 17: Axial velocity cross-sectional pro-files for initial (black line) and op-timized (red line) shape at Z =0 . m . Location of the crosssection is schematically shown atbottom left. U z ( m / s ) r (m) Initial shapeOptimized shape Figure 18: Axial velocity cross-sectional pro-files for initial (black line) and op-timized (red line) shape at Z =0 . m , corresponding to the fi-nal nodes before the sudden ex-pansion of the initial shape. -6 -5 -5 -5 -5 -5 H e m o l ys i s ( H ) r (m)Initial shapeOptimized shape Figure 19: Hemolysis cross-sectional profilesfor initial (black line) and opti-mized (red line) shape at Z =0 . m . Location of the crosssection is schematically shown atbottom left. -6 -5 -5 -5 -5 -5 H e m o l ys i s ( H ) r (m)Initial shapeOptimized shape Figure 20: Hemolysis cross-sectional profilesfor initial (black line) and opti-mized (red line) shape at Z =0 . m , corresponding to the fi-nal nodes before the sudden ex-pansion of the initial shape.21 E-081E-071E-061E-05 GW HO ZT
22 % 19 % 20 %22 % 20 % 20 % H i nde x Set of ParametersInitial shape, k = 1Optimized shape, k = 1Initial shape, k = 2Optimized shape, k = 2
Figure 21: Computed values of H index (Eqn. 9) fordifferent set of parameters (cf. Tab. 1)and different k values (Eqn. (cf. 4)). Foreach set of parameters and k value thepercentage decrease of the objective onthe initial w.r.t to the optimized shape isshown next to points. Y-axis in logscale. Shape Param. k H index Initial GW 1 0.460E-06Optimized GW 1 0.359E-06Initial GW 2 0.106E-05Optimized GW 2 0.826E-06Initial HO 1 0.169E-07Optimized HO 1 0.136E-07Initial HO 2 0.337E-07Optimized HO 2 0.270E-07Initial ZT 1 0.135E-06Optimized ZT 1 0.109E-06Initial ZT 2 0.271E-06Optimized ZT 2 0.215E-06Table 3: Exact values from Fig.21.22
Conclusions and Outlook
This paper discusses a continuous adjoint approach for shape optimization targeting tominimize flow-induced hemolysis in biomedical applications. A detailed derivation ofthe adjoint model which stems from the original power-law formulation for hemolysisprediction is reported.A benchmark problem is examined to verify the numerical implementation of the primaland adjoint hemolysis equations. The numerical results compare favorably against thosededuced from analytical solutions of the problem. The validity of the derived sensitivityderivative is put to the test on a two-dimensional stenosed geometry on which a FD studyis realized. The continuous computed sensitivity from the adjoint method is comparedagainst the locally evaluated results of the FD study and a fitness is found.The complete optimization method based on the derived adjoint equations is appliedon a three-dimensional geometry, specifically designed to include flow peculiarities, re-lated to blood damage in medical devices. An optimized shape, in 50 shape updates,able to reduce the numerically predicted flow-induced hemolysis by 22% is found. Theperformance of the optimized shape, in terms of blood damage, is then tested for differ-ent hemolysis-related parameters. It is found that in all cases, the optimized geometryoutperforms the initial one.Overall, the reported method poses great potential for minimizing flow-induced hemoly-sis, in silico conditions, due to its computational efficiency and relative ease of numericalimplementation on a standard CFD solver. While the method is derived based on a spe-cific hemolysis prediction model on structures with rigid walls, future work will targetthe extension of the method to different prediction models as well as a Fluid StructureInteraction (FSI) formulation, so that wall elasticity is accounted for.
Acknowledgments
The current work is a part of the research training group “Simulation-Based DesignOptimization of Dynamic Systems Under Uncertainties” (SENSUS) funded by the stateof Hamburg under the aegis of the Landesforschungsförderungs-Project LFF-GK11. Se-lected computations were performed with resources provided by the North-German Su-percomputing Alliance (HLRN). This support is gratefully acknowledged by the au-thors. 23 ppendices
Appendix A
Staying true to the notation of Sec. 2.2.1, the derivation of the directional derivatives I k , with k = 1 , ..,
4, is presented herein. Each sub-integral expression is repeated for theconvenience of the reader. Finally, the process to arrive at the adjoint equations (Eqns.(23)-(25)) is also illustrated.For I : I (ˆ p, u i , Ω) = Z Ω ˆ p (cid:16) − ∂u j ∂x j (cid:17) d Ω . Since I doesn’t directly depend on p and H L it follows δ p I · δp = δ H L I · δH L = 0 , (40)while for the variation in the direction of the velocity components, by changing the index j to i and performing integration by parts, we get δ u i I · δu i = − Z Ω ˆ p ∂δu i ∂x i d Ω= − Z Γ ( δu i n i )ˆ pd Γ + Z Ω ( δu i ) ∂ ˆ p∂x i d Ω . (41)For I : I (ˆ u i , u i , Ω) = Z Ω ˆ u i (cid:16) ρu j ∂u i ∂x j (cid:17) d Ω . Since I doesn’t directly depend on p and H L it follows δ p I · δp = δ H L I · δH L = 0 (42)and for its variation on the direction of the velocity components24 u i I · δu i = Z Ω ˆ u i ρ (cid:16) u j ∂δu i ∂x j + δu j ∂u i ∂x j + δu j ∂δu i ∂x j (cid:17) d Ω= Z Ω ˆ u i ρu j ∂δu i ∂x j d Ω + Z Ω ˆ u j ρδu i ∂u j ∂x i d Ω + neglect z }| {
HoT = Z Γ ˆ u i ρ ( u j n j ) δu i d Γ − Z Ω ( δu i ) ρu j ∂ ˆ u i ∂x j d Ω+ Z Ω ˆ u j ρ ( δu i ) ∂u j ∂x i d Ω= Z Γ δu i ˆ u i ρ ( u j n j ) d Γ − Z Ω ( δu i ) (cid:18) u j ρ ∂ ˆ u i ∂x j − ˆ u j ρ ∂u j ∂x i (cid:19) d Ω . (43)The volume integral on the 5th line of Eqn. (43) could also be integrated by parts. Thiswould lead to an exchange of the adjoint and primal velocities on the volume integraland an additional surface integral. However, neglecting this operation could stabilize thenumerical solution of the adjoint system [29]. In this study, the derivation as presentedin Eqn. (43) is solely considered.For I : I (ˆ u i , u i , Ω) = Z Ω ˆ u i (cid:16) − ∂∂x j (2 µS ij ) (cid:17) d Ω . As before δ p I · δp = δ H L I · δH L = 0 (44) δ u i I · δu i = − Z Ω ˆ u i ∂∂x j (cid:20) µ (cid:16) ∂δu i ∂x j + ∂δu j ∂x i (cid:17)(cid:21) d Ω= − Z Γ ˆ u i µ (cid:16) ∂δu i ∂x j + ∂δu j ∂x i (cid:17) n j d Γ | {z } Γ + Z Ω ∂ ˆ u i ∂x j µ (cid:16) ∂δu i ∂x j + ∂δu j ∂x i (cid:17) d Ω | {z } I ∗ , (45)where I ∗ can be expanded further using a second integration by parts as25 ∗ = Z Ω ∂ ˆ u i ∂x j µ ∂δu i ∂x j d Ω + Z Ω ∂ ˆ u i ∂x j µ ∂δu j ∂x i d Ω= Z Γ ( δu i ) µ ∂ ˆ u i ∂x j n j d Γ − Z Ω δu i ∂∂x j (cid:20) µ ∂ ˆ u i ∂x j (cid:21) d Ω+ Z Γ ( δu j ) µ ∂ ˆ u i ∂x j n i d Γ − Z Ω δu j ∂∂x i (cid:20) µ ∂ ˆ u i ∂x j (cid:21) d Ω= Z Γ ( δu i ) µ (cid:20) ∂ ˆ u i ∂x j + ∂ ˆ u j ∂x i (cid:21) n j d Γ | {z } Γ − Z Ω ( δu i ) ∂∂x j (cid:20) µ (cid:16) ∂ ˆ u i ∂x j + ∂ ˆ u j ∂x i (cid:17)(cid:21) d Ω . (46)Γ + Γ are subject to even further expansionΓ + Γ = Z Γ (cid:20) ( δu i ) µ h ∂ ˆ u i ∂x j + ∂ ˆ u j ∂x i i n j − ˆ u i µ h ∂δu i ∂x j + ∂δu j ∂x i i n j (cid:21) d Γ= Z Γ µn j (cid:18) ∂ ˆ u j ∂x i δu i − ∂δu j ∂x i ˆ u i (cid:19) d Γ + Z Γ µn j (cid:18) ∂ ˆ u i ∂x j δu i − ∂δu i ∂x j ˆ u i (cid:19) d Γ . (47)Expanding the first integral on Eqn. (47) leads to Z Γ µn j (cid:18) ∂ ˆ u j ∂x i δu i − ∂δu j ∂x i ˆ u i (cid:19) d Γ = Z Γ µ ∂ ˆ u n ∂x i δu i d Γ − Z Γ µ ∂δu n ∂x i ˆ u i d Γ= µ ˆ u n δu i n i | {z } δu n ∂ Γ − µ ˆ u i n i |{z} ˆ u n δu n ∂ Γ − Z Γ ˆ u n µ ∂δu i δx i d Γ + Z Γ δu n µ ∂ ˆ u i ∂x i d Γ= 0 . (48)The last two integrals of Eqn. (48) vanish because the primal velocity field is asymp-totically divergence-free in the boundary so its variation must be divergence-free tooand because we showed that the adjoint velocity must also be divergence-free. Finally,simplifying Eqn. (47) through Eqn. (48) and inserting the remaining terms in Eqn. (45)we get δ u i I · δu i = Z Γ ( δu i ) µn j ∂ ˆ u i ∂x j d Γ − Z Γ ˆ u i µn j ∂δu i ∂x j d Γ − Z Ω ( δu i ) ∂∂x j (cid:16) µ ˆ S ij (cid:17) d Ω . (49)For I : 26 (ˆ u i , p, Ω) = Z Ω ˆ u i (cid:16) ∂p∂x i (cid:17) d Ω . It follows that δ u i I · δu i = δ H L I · δH L = 0 , (50)while for the variation in the direction of pressure we can write δ p I · δp = Z Ω ˆ u i ∂δp∂x i δ Ω= Z Γ ( δp )ˆ u i n i d Γ − Z Ω ( δp ) ∂ ˆ u i ∂x i d Ω . (51)Having expressed all sub-integrals in the form of surface and volume integrals we cansuperpose all contributions (considering also those from the main body of the paper)into the directional derivatives of the Lagrange functional so that δ p L · δp = Z Γ ( δp )ˆ u i n i d Γ − Z Ω ( δp ) ∂ ˆ u i ∂x i d Ω (52) δ u i L · δu i = Z Γ ( δu i ) (cid:20) − ˆ pn i + ρ ˆ u i u j n j + µ ∂ ˆ u i ∂x j n j + H L ˆ hn i − ˆ hB ij n j + ∂j Γ ∂u i (cid:21) d Γ+ Z Ω ( δu i ) " ˆ u j ρ ∂u j ∂x i − u j ρ ∂ ˆ u i ∂x j − ∂∂x j (cid:16) µ ˆ S ij − ˆ pδ ij (cid:17) + ∂ (ˆ hB ij ) ∂x j − H L ∂ ˆ h∂x i d Ω − Z Γ µ ˆ u i ∂δu i ∂x j n j d Γ (53) δ h L · δh = Z Γ ( δh ) (cid:16) ˆ h ( u j n j ) + ∂j Γ ∂H L (cid:17) d Γ + Z Ω ( δh ) (cid:16) u j ∂ ˆ h∂x j − C β τ ab ˆ h (cid:17) d Ω . (54)Finally, by demanding that the integrals of the form of Eqn. (22) disappear for everytest function, δy , one arrives at the adjoint equations (Eqns. (23) - (25)) as well as theirBCs (Eqns. (26)). It is worth noting that the second surface integral on the RHS ofEqn. (53) does not conform to the integrals of Eqn. (22). This term might contributeto the sensitivity derivative depending on the adjoint BCs.27 ppendix B Additional information for the computation of surface sensitivity
The mathematical manipulations that are followed for the derivation of Eqn. (31) fromEqn. (30) are presented herein. The latter equation is rewritten here for the convenienceof the reader δ c L · δc = − Z Ω ˆ y (cid:16) δ y R · δy (cid:17) d Ω . The RHS of the equation can also be written as