Adjoint-based Shape Optimization for the Minimization of Flow-induced Hemolysis in Biomedical Applications

This paper reports on the derivation and implementation of a shape optimization procedure for the minimization of hemolysis induction in biomedical devices. Hemolysis is a blood damaging phenomenon that may occur in mechanical blood-processing applications where large velocity gradients are found. An increased level of damaged blood can lead to deterioration of the immune system and quality of life. It is, thus, important to minimize flow-induced hemolysis by improving the design of next-generation biomedical machinery. Emphasis is given to the formulation of a continuous adjoint complement to a power-law hemolysis prediction model dedicated to efficiently identifying the shape sensitivity to hemolysis. The computational approach is verified against the analytical solutions of a benchmark problem and computed sensitivity derivatives are validated by a finite differences study on a generic 2D stenosed geometry. The application included addresses a 3D ducted geometry which features typical characteristics of biomedical devices. An optimized shape, leading to a potential improvement in hemolysis induction up to 22%, is identified. It is shown, that the improvement persists for different, literature-reported hemolysis-evaluation parameters.


Introduction
The ever-growing advances in medicine, engineering and material science have led to the development of biomedical devices, such as blood pumps, which allow long-term patient care and significantly improve quality of life. Despite all the advances, a critical task for the design and development of such devices [33,38] is still the minimization of shear-induced blood damage (i.e. hemolysis) to guarantee good bio-compatibility.
In the context of this manuscript, hemolysis refers to the mechanical damage of red blood cells due to excessively high stress induced by peculiarities of the blood flow. It can 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 immune system [4]. Hemolysis induction is encountered in many biomedical devices, where large velocity gradients are found [3], as well as in vivo conditions when vessels delivering blood are kinked or stenosed [9]. The shape of the respective artificial devices or vessels is believed to play a crucial part in the induction of blood damage due to its decisive fluid dynamic role. A computationally efficient shape optimization framework would be therefore appreciated for the development of next-generation biomedical machinery.
The success of a numerical optimization process partially depends on the accuracy of the hemolysis-prediction model. While the study of shear-induced hemolysis has been of interest for many experimental, in vitro conditions for quite some time [6,12,39], it is also becoming increasingly important in silico conditions [30,37]. Computational fluid dynamics (CFD) offer the possibility to predict blood damage in a purely numerical manner and employ the damage prediction model as a cost function in an optimization framework. A variety of such damage prediction models have been proposed [37]. Most of them relate hemolysis with the magnitude of the shear stress and an exposure period using 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 blood damage in a variety of flows has not yet been established [9]. The present manuscript does not aim to advocate the merits of a specific hemolysis evaluation or hemolysis prediction model, but is mainly concerned with the adaptation and integration of a classical model into an adjoint shape optimization framework. Employing an Eulerian hemolysis-prediction model alongside a common CFD solver, we formulate a Partial-Differential-Equation (PDE) constrained optimization problem, which targets to minimize 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-based optimization, multiple ways, ranging from stochastic [1,25] to deterministic [2,34] optimization methods, could be followed. The present research is concerned with the efficient computation of the derivative of the objective with respect to (w.r.t) the control. The aforementioned derivative is subsequently used by a deterministic gradient-based steepest descent method, which drives the controlled shape towards an improved state. To that 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] since the pioneering works of Pironneau [26] and Jameson [15], due to its superior computational efficiency. In specific, the attractiveness of the method lies on the fact that the computation of the derivative is independent from the size of the control. For further details on the merits and drawbacks of the adjoint approach, the interested reader is referred to [7,27].
The remainder of this paper is organised as follows: Section 2 presents the mathematical model of the primal and adjoint problem. The same section also reports on the employed numerical procedure and utilized boundary conditions (BCs). In Section 3, a benchmark case is studied to verify the code reliability. Moreover, the accuracy of the computed 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 a full shape optimization process are presented and their dependency on parameters of the nonlinear hemolysis evaluation model is discussed. The paper closes with conclusions and outlines further research in Section 5.
Within this publication, Einstein's summation convention is used for repeated lower-case Latin subscripts. Vectors and tensors are defined with reference to Cartesian coordinates.

Mathematical Model
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 equations is discussed. A presentation of the newly developed adjoint model then follows with detailed discussions on adjoint-hemolysis specific points of interest.

Primal Flow
Throughout this paper, blood is treated as an incompressible, Newtonian fluid. The assumption of non-Newtonian behaviour has been shown to satisfactorily predict the flow of blood in tubes with diameter ranging from 130 to 1000 µm [24]. However, all the applications considered herein refer to ducted systems with a larger diameter and thus the Newtonian assumption is preferred. Furthermore, we assume all applications in laminar and steady-state conditions. The blood flow in a domain (Ω) follows from the Navier-Stokes (NS) equations for the conservation of volume and momentum, viz.
Here u i , p, S ij = 0.5(∂u i /∂x j + ∂u j /∂x i ), δ ij , ρ and µ refer to the components of the fluid velocity vector, static pressure, components of the strain-rate tensor, Kronecker delta 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 a prediction or evaluation model in a one-way coupling with the CFD solver [37]. This implies that the hemolysis model receives information from the NS equations (1,2) with no retro-action on the employed fluid properties (ρ, µ). The hemolysis prediction model considered in this study originates from the power-law equation, first introduced by Giersiepen et al. [6], and reads The hemolysis index H denotes a measure for the released hemoglobin to the total hemoglobin 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 is exposed to the stress. The remaining constants (C, α, β) are introduced to fit experimental data. (3). The notation corresponds to the initials of the investigators that performed the experiments (GW [6], HO [12], ZT [39]). The final two columns correspond to the type of blood that was used and to the maximum stress that was applied during the experiments.
invariant of the stress tensor τ ij = 2µS ij is frequently used to determine the scalar stress representation τ in combination with a user specific parameter k = 1, 2, 3 Due to tr(τ ij ) = tr(2µS ij ) ∼ ∂u i /∂x i , the first contribution to I τ 2 in (4) vanishes for incompressible Newtonian fluids, which yields To incorporate Eqn. (3) into an adjoint optimization framework, a PDE-based hemolysis prediction is advantageous. To this end, a linearization strategy is used [5] for (3), i.e.
Substituting H L = H 1 β , the sought material derivative expression follows from (6) Assuming steady-state, a modification of the residual form of (7) reads Note that (8) is enhanced by (1 − H L ) to ensure H L < 1. Equation (8) will be referred to as (primal) hemolysis equation throughout this manuscript and serves to formulate the objective functional of the optimization. The merits of expression (8) refer to its straightforward derivation from the original power-law formulation (3) and its suitability for showcasing the adjoint method in the context of blood damage. While it is true that hemolysis does not occur for small values of τ < τ th , where τ th corresponds to some threshold value (which might also be related to the exposure time), this aspect is deliberately ignored in the present effort. Furthermore, the footing of the hemolysis model assumes a homogeneous stress representation. Even though this is a strong assumption for most practical applications, the predictions have been shown to be sufficiently reliable when relatively comparing different geometries [32]. In the present context of shape optimization, it is thus fair to assume that the qualitative information on a reduction of the induced hemolysis is reliable.
Equation (8) 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-functional in this manuscript.

Adjoint Method
This subsection outlines the formulation of the adjoint problem at hand which is initially expressed as an optimization problem constrained by a set of PDE, namely R(y, c) = 0, viz. min Here J is the objective or cost functional, y is the vector of primal state variables, which consists of p, u i and H L , and c is the control parameter which we need to find in the set of admissible states C ad . In this study, c refers to the shape of the structure in which blood flows. It is thus a subset of ∂Ω or Γ -the latter notation is preferred throughout this manuscript-and can be further restricted by considering only specific sections of the structure, namely Γ D ⊂ Γ, where the index D denotes design.
Based on (9) the objective functional under investigation can be written as while the set of PDE, R(y, c) = 0, consists of Eqns. (1), (2) and (8). Having formulated the optimization problem as a constrained problem in Eqn. (10), the Lagrange principle to eliminate the constraints by employing appropriate Lagrange multipliers is utilized [29]. In this context, a continuous Lagrange functional is defined as are required to be zero, it is apparent that L is equal to J. The optimization problem (10) is, thus, equivalent to minL(ŷ, y, c), whereŷ is the vector of adjoint variables and y is no longer constrained. For a minimum of (12) the total variation δL vanishes The terms FD1, FD2 and FD3 denote the variation of the Lagrangian in the direction of adjoint (δŷ) and primal (δy) state as well as the control (δc), respectively. They can be computed by utilising the functional derivative into the respective direction. FD1 leads to the known set of primal equations that need to be satisfied for every control state. FD2 yields the accompanying adjoint equations. Finally, FD3 gives rise to a sensitivity 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 following three optimality conditions are obtained for vanishing FD2 A detailed overview on PDE constraint optimization can be found in [13].

Adjoint Flow Equations
Each optimality condition (14) leads to a PDE that needs to be satisfied in order for FD2 to vanish. The essence of the adjoint method is to identify an adjoint stateŷ that satisfies (14), so that the total variation δL depends only on the variation of the control (δ c L · δc). For the purpose of this paper it is deemed beneficial to split (12) into sub-integrals consisting of isolated terms of the primal equations, viz.
The computation of the functional derivatives (14) requires the computation of the derivatives of each sub-integral I k , with k = 1, 2..., 6, as well as the derivatives of the objective functional (11). The latter involves δ p J · δp = 0 and Note that the objective functional and its directional derivatives only exist at the outlet of the domain, which is not a design surface (Γ D ). The contribution of the objective functional 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 inlet and the outlet pressure, the dual flow is driven by the relative difference of the hemolysis index.
The calculation of the variations δ y I k · δy, with k = 1, 2.., 4 has been previously shown in many manuscripts [22,29]. It is therefore omitted from the main body of this paper and can be found in App. A. The derivation for the hemolysis relevant terms I 5 and I 6 represents a core contribution of this paper. For δ y I 5 · δy, we can utilize Gauss's divergence theorem and obtain Although Eqn. (17) doesn't necessarily require integration by parts to be included in the adjoint equations, the respective volume integral is expanded in the context of this work for computational purposes.
Following the same strategy for I 6 , yields δ p I 6 · δp = 0. As for δ u i I 6 · δu i , we first define as a scalar quantity not subjected to any variation w.r.t the velocity. Using the scalar stress expression (5), the integral I 6 is expressed as and its linearized variation in the velocity direction reads where (19) is further expanded via integration by parts as Finally, δ H L I 6 · δH L is calculated as Having expressed the variation of J and I k in terms of several boundary and volume integrals, we can superpose the variation of the Lagrange functional as the sum of one boundary and volume integral, viz.
The adjoint equations as well as their corresponding boundary conditions arise by demanding that the integrands in Eqn. (22) vanish for any test function δy. Therefore, following the aforementioned methodology for the optimality conditions (14), the resulting constraints in the interior of the domain (Ω) read Equations (23) as well as (24) represent the adjoint companions to the continuity and momentum equation, respectively. In this sense, Eqn. (25) is referred to as adjoint hemolysis equation throughout this manuscript. The adjoint momentum equation is enhanced with the last two terms on the left-hand-side (LHS) of the equation that include contributions from adjoint and primal hemolysis. This is due to the fact that on the primal side of the system, the velocity appears twice in the hemolysis equation, once in the advection and once in the form of a gradient in the source term. Both terms act 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, the adjoint hemolysis equation doesn't require any information from the solution of the adjoint continuity and momentum equations. Once again, we have a one-way coupling but this time the direction of the coupling and thereby the algorithmic order of sequence is reversed. Practically this means that the adjoint hemolysis equation could be solved after the solution of the primal equations and prior to the solution of the other two adjoint Eqns. (23) & (24).
The adjoint system (23)(24)(25) is an extension of a classical adjoint (incompressible) NS system, which corresponds to (23)(24) excluding the two hemolysis terms. This is similar to the adjoint complement of the Volume of Fluid (VoF) method [16,18], used for multiphase applications, and enables a straightforward implementation, provided that an adjoint multi-phase solver exists. Although the adjoint equations are per definition linear, several cross coupling terms might introduce a severe stiffness to the densely coupled PDE system. A viable workaround to stabilize the iterative procedure refers to the introduction 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

Computation of Surface Sensitivity
Ensuring a vanishing residual of the primal and adjoint system of equations leads to vanishing FD1 and FD2 terms of (13). The FD3 term allows for the computation of a surface sensitivity which is used by a gradient-based algorithm to drive the shape towards an improved state. Following (12), FD3 is expressed as Since the primal flow equations (1,2,8) are also fulfilled under the total variation, one obtains δR(y, c) = δ c R · δc + δ y R · δy = 0.
Therefore, the variation w.r.t. the control transforms to Furthermore, based on the formulation of our objective functional J, which is defined only on the outlet, δ c J ·δc = 0 since c ≡ Γ D and Γ D ∩Γ out = ∅. The sensitivity derivative can thus be written as The right-hand-side (RHS) of (30) is developed further and by taking into consideration the BCs of the problem as well as a Taylor expansion of the velocity w.r.t a perturbation of the design wall (see App. B or [20]), the sensitivity to be computed reads where u i(t) ,û i(t) denote the part of the primal and adjoint velocity tangential to the surface, respectively and n denotes the normal to the surface. Interestingly, the shape derivative 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), and thereby propagate into the shape derivative through the adjoint velocity gradient.

Boundary Conditions
In what concerns the application studied in this work, the flow boundary Γ consists of three parts, namely the inlet (Γ in ), outlet (Γ out ) and wall (Γ W ). We split the wall boundary into two portions, Γ D and Γ B = Γ W \Γ D , where the first one involves the parts under design while the latter one is bounded to the initial configuration. The complete boundary domain is thus described as While the BCs of the primal system of equations are selected based on the physical properties of the problem, the adjoint BCs are deduced based on the necessity of vanishing the surface integrals of Eqn. (22). In many boundary patches this necessity is inherently satisfied due to the primal BCs, though in general it reduces to satisfying Eqns. (26). Furthermore, it is worth mentioning that the BCs strongly depend on the objective functional, if that is defined in a subset of Γ. The complete set of BCs for the primal and adjoint problem are summarized in Tab. 2.
Boundary Patch

Numerical Procedure
The numerical procedure for the solution of the primal and adjoint system is based upon the Finite Volume Method (FVM) employed by FreSCo + [28]. Analogue to the use of integration-by-parts in deriving the continuous adjoint equations [19,21], summation-byparts is employed to derive the building blocks of the discrete adjoint expressions. A detailed derivation of this hybrid adjoint approach can be found in [16,18,31]. The last two terms of the adjoint momentum equation, involving hemolysis contributions, are added explicitly to the RHS. The segregated algorithm uses a cell-centered, collocated storage arrangement for all transport properties. The implicit numerical approximation is second order accurate in space and supports polyhedral cells. Both, the primal and adjoint pressure-velocity coupling is based on a SIMPLE method and possible parallelization is realized by means of a domain decomposition approach [35,36].

Optimization Framework
The complete shape-optimization framework is summarized in Fig. 1. The efficiency of the method is illustrated by taking into consideration the required computational budget for one complete optimization process based on the presented flowchart. Assuming that the cost for solving the primal problem is equal to 1 equivalent flow solution (EFS) then the cost for solving the adjoint one is approximately also equal to 1 EFS. The final two processes (represented by rectangles on the flowchart) have a practically negligible cost compared to the EFS unit. The complete optimization cost is thus equal to i max times 2 EFS, regardless of the size of the control which for our case is proportional to the number of discretized nodes on Γ D . In the context of this study, the steepest descent method is considered to advance the shape.
The computed surface sensitivity S L might possibly be rough. We thus employ the Laplace-Beltrami [17] metric to extract a smooth gradient G L i through a numerical approximation of where ∆ Γ refers to the Laplace-Beltrami operator (∆ Γ = ∆ − ∆ n ), λ corresponds to a user-defined control of smoothing and n i denotes the normal vector at each face of Γ. Update shape to c i+1 using the steepest descent method.
Is i < i max ? End YES NO Figure 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 from Eqn. (33), where q = 1/(W D + ǫ), W D is the wall normal distance and ǫ = 10 −20 .

Verification and Validation Studies
This section verifies the implementation of the primal and adjoint hemolysis system for 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.

Verification
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 skipped and the reader is referred to [10]. The solution of the primal flow reads where r corresponds to the radial direction, R marks the pipe radius, U max refers to the centerline velocity, z to the axial coordinate and the direction of the primary flow velocity u z .
The considered verification case is sketched in Fig. 2. The entrance and exit planes of the pipe refer to z = 0 and z max = 2 m. The pipe radius reads R = 0.5 m. The fully developed axial velocity profile of the laminar flow utilizes U max = 10 m/s. The fluid properties are selected based on a prescribed Reynolds number of Re = 2U R/ν = 2000, where U and ν refer to the bulk velocity and kinematic viscosity, respectively. The three-dimensional geometry is discretized with approximately 75k control volumes on a 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 solution can be stated using the analytical primal flow solution (34), viz.
If we abbreviate Λ = C 1 β τ α β , then the solution of (35) readŝ The integration constant K is calculated based on the BCs of the adjoint hemolysis equation which results from satisfying the final expression in (26) and demandsĥ| out = −ρβH β−1 L | out . Finally, the adjoint hemolysis field in a fully developed pipe flow readŝ 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 values as calculated by Eqn. (34). The comparison is realized for all three sets of hemolysis model parameters outlined in Tab. 1. As can be seen, all computed values are fitting the analytical solutions to a satisfying degree. Figure 5 compares the computed adjoint  We would like to remark on the nature of the adjoint hemolysis equation. As can be seen in Fig. 5, the adjoint hemolysis profile changes only slightly for different axial positions. The BC at the outlet, which readsĥ| out = −ρβH β−1 L | out , dominates the complete field since all three cases correspond to β < 1 and H L << 1. To avoid numerical errors that would arise for H L | out = 0, the BCs on the outlet is reformulated toĥ| out = −ρβ(H L + ǫ) (β−1) | out , where ǫ = 10 −20 . Based on the previous comments, one could assume that applying the bulk outlet adjoint hemolysis profile to the whole field suffices. 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 description of the complete model.

Sensitivity Validation
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 accuracy of the computed sensitivity. In order to realize this, the computed shape derivative is compared against locally evaluated second order accurate finite differences [14,17], viz.
Here c i represent discrete points of the control, ǫ is the magnitude of the perturbation and n i is the normal vector at c i . In practise, the study is realized by deforming the boundary faces of the discretized geometry into their normal direction with a magnitude equal to ǫ. The local boundary perturbations are then transported into the domain based on (33). Figure 6 presents a schematic of the considered symmetric geometry. In specific, the design section (Γ D ) follows from y = A sin π L x − π 1000 4 , where A is the height of the bump, L is the length of the design surface and L/A = 20. The domain is discretized with approximately 70k control volumes on a structured grid which is refined in normal as well as tangential direction towards the wall and Γ D . A detail of the utilized grid 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.3 m/s. The fluid properties (ρ, µ) are assigned to unity to simplify the study. The results of the FD study are shown in Fig. 8 for two perturbation magnitudes, ǫ, at 20 selected discrete positions. The overall agreement between the computed adjoint shape derivative and the FD results is very good. Furthermore, as shown in Fig. 9 the computed objective functional, on the perturbed shapes, exhibits a linear behaviour w.r.t the perturbation size.
Overall, the study validates, on a preliminary basis, the accuracy of the adjoint method, as presented within this manuscript.

Application
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 a 3D geometry. The latter corresponds to a benchmark model, designed by a technical committee to include flow phenomena related to blood damage in medical devices [30].

Investigated Configuration
The geometry under consideration shares characteristics of many blood-carrying medical devices. It includes sections where the flow is accelerating or decelerating, where recirculation and significant variations in shear stress occur. All of these phenomena are believed to be related to blood damage in medical devices [30]. A merit of the adjoint method is that it stems from the primal (physical) problem and does not operate as a black box. It is therefore capable of shedding light into areas of potentially elevated blood-damage capacity. Areas on which an accumulated shape sensitivity is identified are 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 again split into 2 parts, Γ D and Γ B , as was done in the previous section. The inlet and outlet tubes are considered to be bounded, since their shape is believed to be trivial to the blood damage capacity. The remaining structure is classified as Γ D and presented in blue.   [30]. Presented in blue are the sections free for design. The remaining structure is bounded to its initial configuration. All dimensions in (mm). Bottom: Detail of the computational grid on a longitudinal section of the geometry near the design surface.
The geometry is discretized with approximately 1 million control volumes on a butterflylike 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, the grid 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 3 and 3.5 cP , respec-tively, representing blood under physiological conditions. A fully developed laminar axial velocity profile is prescribed at the inlet of the geometry, based on (34), so that the Reynolds number at the throat reads Re T = 500. As regards the hemolysis-prediction model, 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 terms in the adjoint and primal momentum equations are discretized using the second order accurate central difference scheme while the convective terms are discretized through the higher order Quadratic Upstream [Downstream] Interpolation of Convective Kinematics (QU(D)ICK) scheme.

Shape Optimization Study
The optimization study was performed for 50 design candidates. The employed smoothing parameter value (32) reads λ = d T , where d T is the throat diameter. To avoid large deformations of the geometry and the grid, the local displacement field d i computed from Eqn. (33) was scaled to a user-defined constant maximum valued max , viz.
For the study at handd max /d T = 0.25 · 10 −2 . Figure 13 shows the displacement magnitude on the design surfaces of the structure at hand, after the solution of the first primal and adjoint problem. As can be seen, the displacement is accumulated, almost entirely, at the sudden expansion of the geometry, where the highest values of hemolysis occur. Furthermore, re-circulation is occurring with relatively significant values of upstream mass fluxes after the expansion. This is believed to be one of the causes of hemolysis. Therefore, even though there is no direct contribution from the primal or adjont hemolysis to the shape sensitivity, the necessary information from both fields are preserved through the adjoint velocity, which is directly influenced by both H L andĥ. The optimization history is shown in Fig. 15. It can be seen that the optimization starts converging after approximately 40 shapes. The final shape results in a relative reduction of the objective functional by 22%. A comparison of the initial and the optimized shape is displayed in Fig. 14. The optimization algorithm, proceeded into widening and relatively smoothing the sudden expansion of the geometry. This results in a smaller maximum value of the axial velocity as can be seen in Fig. 16. At the same time, nonetheless, the flow inside the bounded portions of the structure (z ∈ [0 : 0.12] m) remains relatively unchanged as shown on the same figure. This is important in terms of biomedical applications due to their ultimate goal of realizing a sensitive task. However, due to the smoothing of the expansion zone, the velocity profiles are changed in the perturbed section of the geometry (cf. Figs. 17 & 18). The velocity profile of the optimized shape 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 the structure, which is also the area on which the maximum values are identified in the initial shape. A direct comparison of the hemolysis profiles at two cross sections of the   initial and final geometry is presented in Figs. 19 & 20. As described in the previous section, the optimization utilized a specific set of parameters for the primal hemolysis model. It is interesting, therefore, to examine the performance of the 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 outperforms the initial one in terms of flow-induced hemolysis. In specific, the most significant improvement 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 shows a robustness of the method in terms of user-defined values.

Conclusions and Outlook
This paper discusses a continuous adjoint approach for shape optimization targeting to minimize flow-induced hemolysis in biomedical applications. A detailed derivation of the adjoint model which stems from the original power-law formulation for hemolysis prediction is reported.
A benchmark problem is examined to verify the numerical implementation of the primal and adjoint hemolysis equations. The numerical results compare favorably against those deduced from analytical solutions of the problem. The validity of the derived sensitivity derivative is put to the test on a two-dimensional stenosed geometry on which a FD study is realized. The continuous computed sensitivity from the adjoint method is compared against the locally evaluated results of the FD study and a fitness is found.
The complete optimization method based on the derived adjoint equations is applied on a three-dimensional geometry, specifically designed to include flow peculiarities, related 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. The performance of the optimized shape, in terms of blood damage, is then tested for different hemolysis-related parameters. It is found that in all cases, the optimized geometry outperforms the initial one.
Overall, the reported method poses great potential for minimizing flow-induced hemolysis, in silico conditions, due to its computational efficiency and relative ease of numerical implementation on a standard CFD solver. While the method is derived based on a specific hemolysis prediction model on structures with rigid walls, future work will target the extension of the method to different prediction models as well as a Fluid Structure Interaction (FSI) formulation, so that wall elasticity is accounted for.

Appendices
The volume integral on the 5th line of Eqn. (43) could also be integrated by parts. This would lead to an exchange of the adjoint and primal velocities on the volume integral and an additional surface integral. However, neglecting this operation could stabilize the numerical solution of the adjoint system [29]. In this study, the derivation as presented in Eqn. (43) is solely considered.
For I 3 : As before δ p I 3 · δp = δ H L I 3 · δH L = 0 (44) where I * 3 can be expanded further using a second integration by parts as Γ 1 + Γ 2 are subject to even further expansion Expanding the first integral on Eqn. (47) leads to The last two integrals of Eqn. (48) vanish because the primal velocity field is asymptotically divergence-free in the boundary so its variation must be divergence-free too and 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 For I 4 : It follows that while for the variation in the direction of pressure we can write Having expressed all sub-integrals in the form of surface and volume integrals we can superpose all contributions (considering also those from the main body of the paper) into the directional derivatives of the Lagrange functional so that Finally, by demanding that the integrals of the form of Eqn. (22) disappear for every test function, δy, one arrives at the adjoint equations (Eqns. (23) -(25)) as well as their BCs (Eqns. (26)). It is worth noting that the second surface integral on the RHS of Eqn. (53) does not conform to the integrals of Eqn. (22). This term might contribute to the sensitivity derivative depending on the adjoint BCs.
Substituting Eqn. (58) to (57) we get where the last term of Eqn. (59) vanishes because the primal continuity equation asymptotically holds on the boundaries of the domain also. Finally, if we consider the velocity components as the summation of their normal and tangential parts, one arrives at the expression of Eqn. (31).