Adjoint shape sensitivities of blood flows considering non‐Newtonian properties

This article discusses the derivation and numerical implementation of an adjoint system, to the primal Navier–Stokes equations, for the computation of shape sensitivities of ducted blood flows considering non‐Newtonian fluid properties. The ever‐growing advancements in blood flow simulations are, naturally, accompanied by an increased interest in the optimization of related medical devices. In the majority of the computational studies, the Newtonian assumption is used to describe the rheology of blood. While this assumption has been shown to satisfactorily capture the flow when it is governed by high shear rates, it falls short at low shear rates. A rich variety of viscosity models has been proposed to tackle this shortcoming. In this article we show how such models can be incorporated into an adjoint system targeting to produce the shape sensitivity which can be used by a gradient‐based optimization method for the minimization of an objective functional. A general formulation of the adjoint equations is proposed, in which contributions of the non‐Newtonian properties explicitly occur. The numerical implementation is discussed and the validity of the method is assessed by means of numerical experiments of steady blood flows in a 2D stenosed duct, where results are compared against second‐order finite‐difference (FD) studies. The proposed methodology is then applied to CAD‐free, gradient‐based shape optimizations of an idealized 3D arterial bypass‐graft operating at three relevant Reynolds numbers. It is observed that the impact of the adjoint viscosity treatment is amplified in low shear‐rate flow regimes while fades for higher shear‐rates, analogous to its primal counterpart.


INTRODUCTION
Computational fluid dynamics (CFD) has become an indispensable tool for many industries and fields, concerned with blood flows.On the one hand, it has been extensively used for the development of medical devices, such as ventricular assist devices, 1-3 stents, 4 hemodialysis cannulas 5 and more.The extent to which CFD has positively permeated the medical device industry is highlighted by the interlaboratory study, conducted by the U.S. Food and Drug Administration, for the assessment of CFD simulations on an idealized medical device. 6On the other hand, it has also served to further our fundamental understanding of circulatory diseases and effective treatment of such. 7For example, predictions of the resulting hemodynamics after an invasive surgical procedure for the implantation of prosthetic components are crucial to the success of the treatment.However, such predictions are inevitably difficult to make in vitro, that is, in the laboratory, or in vivo, that is, on a living patient, conditions.To this end, CFD has enjoyed great appeal, being able to overcome restrictions arising in the real world while also offering credible estimations for complex hemodynamics, see for example, Reference 8 and references therein.Overall, computational techniques have been widely used, in recent years, for a vast variety of applications involving blood flow.While each application poses its own difficulties and complexities, a common challenge for all is the appropriate modeling of the rheological behavior of blood.
Strictly speaking, blood is a non-Newtonian fluid due to the deformable red blood cells in plasma and its overall complex nature.However, more often than not, the Newtonian assumption is employed to simplify the complexity of the investigated problems.This assumption has been shown to satisfactorily capture the rheology of blood in flows where shear rates are high, while falling short in flows of low shear rates due to the shear-thinning behavior of blood. 9Several non-Newtonian models have been proposed to deal with this inadequacy and account for pseudo-plastic behavior.Some of such models are the Power-law, 10 Casson's law 11,12 -and modified versions of it 13 -and Carreau's law. 14Even though these models have frequently been used for the numerical investigation of blood flows in a variety of geometries, they still cannot describe blood in a generalized context as they cannot account for viscoelastic properties.Nevertheless, many computational studies of blood flow have employed them due to their straightforward description as well as implementation in a numerical framework.
Concurrently, the ever-growing advances in computational sciences have not only enhanced the studies of blood flows but also enabled the optimization of biomedically relevant shapes.In this context, one searches for optimal configurations of the domain so that a cost or objective functional is minimized.The aforementioned procedure can be realized with various methods, ranging from stochastic 15 to deterministic, 16 with the latter usually requiring the gradient of the objective with respect to (w.r.t.) the shape.Depending on the approach for obtaining this gradient, the computational cost can be restrictive, in applications of complex 3D shapes, for example, finite differences (FD) method, or proportional to the cost of a flow solution, for example, adjoint method. 17Due to its computational efficiency, therefore, the use of the adjoint method has expanded from traditional CFD-based optimization problems, such as optimization of transportation vehicles, [18][19][20] to applications of biomedical engineering. 21,22However, in the majority of biomedical cases, in which adjoint-based shape optimization is considered, blood is modeled exclusively as a Newtonian fluid in either primal and adjoint or viscosity derivatives are dropped from the adjoint equations. 23The latter is similar to the simplifying assumption of frozen viscosity that one might find in the adjoint (dual) problem of turbulent flows, for example, see References 24-27.This article is concerned with the derivation of a consistent dual problem to the Navier-Stokes (NS) equations with a generalized Newtonian constitutive model.Adjoint-based shape sensitivities for the optimization of ducted blood flows are proposed and implementation aspects are discussed.
The remainder of this article is organized as follows: In Section 2, we present the derivation of the adjoint model and highlight the differences w.r.t. the traditional adjoint equations of a Newtonian fluid.The problem is further specified for internal flow configurations and two objective functionals.Section 2 closes with a formulation of the shape sensitivity and a discussion on the employed numerical procedure.Section 3 assesses the accuracy of the shape sensitivity in the context of a FD study.Therein, we also discuss the impact of neglecting non-Newtonian contributions on the adjoint problem.Subsequently, in Section 4, the adjoint model is employed to enable gradient-based shape optimization studies of a 3D idealized arterial bypass-graft.Comparisons are made between the complete non-Newtonian adjoint model and a reduced one, corresponding to the "frozen" viscosity assumption, for three Reynolds numbers.The optimized shapes are then assessed through unsteady simulations to evaluate their efficiency in more realistic conditions.The article closes with conclusions and outlines future directions 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.

METHODOLOGY
This section is dedicated to the formulation and numerical implementation of the adjoint problem for an incompressible fluid in laminar, steady-state flows.We consider three prominent non-Newtonian fluid models, frequently used to describe the shear-thinning behavior of blood, that is, the Power-law (PL), the modified-Casson (MC) and the Carreau models (C).Each of these models describes the fluid viscosity by a non-linear relation to the second invariant of the strain-rate tensor.
In this work the shear-rate ̇, is defined as (1) Based on the above, we opt to define the effective viscosity as μ = F( ̇).Table 1 summarizes the considered effective viscosity expressions for the study at hand.In the context of this work, we employ the constants for each model as presented in Table 1.These correspond to frequently used values for the modeling of the shear-thinning behavior of blood, see for example, References 10 and 13.These values are assumed to be independent of the flow and thus constant throughout a simulation.As regards the (PL) model, we also bound the effective viscosity by μ ∈ [0.002, 0.05] Pa ⋅ s.

Adjoint method
The steady-flow domain is denoted by Ω while its boundary refers to Γ.In what concerns the applications studied in this work, Γ consists of three parts, namely the inlet (Γ in ), outlet (Γ out ) and wall (Γ W ). In the context of shape optimization, the latter is usually split into a section subjected to design (Γ D ) and a section bound to its initial configuration (Γ B ).The complete boundary domain is described as Let J(y, c) denote the objective functional to be minimized, which depends on the fluid dynamic state variables y and the control parameter c.The optimization problem can be stated as min c∈C ad J(y, c) s.t.R(y, c) = 0. (2) In our case, y = [p, u i , μ], where p is the pressure and u i marks the cartesian velocity components.The control parameter c refers to the shape subjected to the optimization and is thus directly related to Γ D .The constraints appearing in (2), R = [R p , R u i , R μ ] T , refer to the residual form of the governing equations that must hold in Ω, namely, TA B L E 1 Constitutive models for the description of effective viscosity.

F( γ)
List of constants Units  Here  ij and  refer to the Kronecker delta components and the fluid's density, while F( ̇) is obtained from Table 1 for the respective viscosity model.The constraints of the considered problem are (a) the NS equations as well as (b) the closure equation of the non-Newtonian viscosity.We therefore introduce a continuous Lagrange functional as where ŷ = [p, û i , μ] are the Lagrange multipliers.
The following optimality conditions are used to build the adjoint system, namely, As a starting point for the derivation, a general description of J that distinguishes between boundary contributions along Γ and volume specific, interior contributions in Ω, as proposed in Reference 26, is given If one splits the volume integral in Equation ( 6) into subintegrals consisting of isolated terms of the primal equations, namely, ) dΩ

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟
, where in accordance with Equation ( 6) the computation of the optimality conditions (7) can be written as The expressions  y I k ⋅ y for y ≡ p, u i and k = 1, … , 4 have previously been derived in many papers, see for example, Reference 26, and details are omitted here for reasons of brevity.Furthermore, it is evident that  μ I k ⋅  μ identically vanishes for k = 1, 2 and 4. The main text is therefore concerned with the derivation and nature of  μ I 3 ⋅  μ and  y I 5 ⋅ y which correspond to the core contributions of this article.By utilizing Gauss's divergence theorem one obtains where n j denotes the normal to the surface vector.For the variations of I 5 one obtains: TA B L E 2 Tensor X ij for each viscosity model.The units of its components are Here, X ij is a symmetric tensor proportional to the strain-rate tensor S ij .The derivation of Equation ( 12) is outlined in Appendix A and summarized in Table 2 for each investigated non-Newtonian model.The optimality conditions (7) can now be written as Here, Ŝij = 0.5(û i ∕x j + û j ∕x i ) refers to the adjoint strain-rate tensor.Boxed terms denote additional contributions w.r.t. the classic optimality conditions for a Newtonian fluid.

Adjoint equations
The adjoint equations arise by demanding that the integrands of the volume integrals of Equations ( 14)-( 16) vanish for any y, namely, Boxed terms denote additional contributions w.r.t. the classic adjoint NS equations for a Newtonian fluid.In line with its primal companion, Equation ( 19) is algebraic (here w.r.t.adjoint viscosity) and therefore one can write where Γij denotes the antisymmetric part of the adjoint velocity gradient and thus its inner product with the symmetric S ij vanishes.
At this point we would like to remark on the nature of the adjoint equations ( Equations 17,18,19/20).First of all, one can notice that for a Newtonian fluid, in which X ij = 0 (cf.Table 2), Equations (17) and (18) lead to the classical adjoint (incompressible) NS system for a generalized objective as described by Equation (8).Furthermore, Equation ( 16) is inherently satisfied since μ = C (cf.Table 1) and  μ = 0 both in the flow domain and its boundaries.Therefore, an adjoint viscosity equation is obsolete.Due to the similarities between the generalized primal and adjoint equations, the computational implementation of the latter can follow a similar algorithmic strategy to that of the primal equations.
Finally, we would like to state that due to the algebraic nature of Equation ( 5) one may arrive at the same adjoint equations by substituting μ in Equation (4) from Equation (5) and omitting I 5 in Equation (6).However, the derivation then becomes significantly more cumbersome and is therefore not shown herein.

Objective functional
This work considers two objective functionals to showcase the generalized adjoint NS system.First, we focus on the minimization of dissipated power, due to the popularity of this objective in optimization studies of ducted flow systems. 26,28e then consider an objective that directly relates to the effective viscosity field and to shear stresses.The latter is studied due to its clinical relevance in applications of blood flow.A further distinction between the two is made based on their habitat, that is, where they are defined.In specific, the former objective is defined on the fluid domain's boundary (Γ) while the latter on an observation region of the domain (Ω o ⊂ Ω), thus resulting to contributions of different nature in the adjoint system.The exemplary application presented in this work refers to the shape optimization of an arterial bypass graft.For such cases, alternative objective functionals of biomedical interest are for example, the oscillatory shear index 29 or the spatial minimum or maximum wall shear stresses over a cardiac cycle. 30However, these metrics require transient optimizations, which go beyond the scope of this article.

Objective 1: Dissipated power
The motivation for considering the minimization of dissipated power as an objective is twofold.Firstly, it refers to a frequently used objective functional in adjoint-based optimization.It is therefore showcased herein for the validation of the method (see Section 3) and supports the reproducibility of the presented method with minimum alterations to existing adjoint codes.Secondly, its biomedical impact is not trivial.Studies have suggested the use of energy dissipation as an index to quantify the severity of aortic stenosis 31 or discussed the metric as a crucial parameter in certain surgical techniques, 32 such as Fontan-type interventions. 33he computation of dissipated power follows from the net inward flux of total pressure through the fluid domain's boundaries, namely, It then holds, The contributions of this objective thus enter the adjoint system exclusively through the boundary integrals of Equations ( 14)- (16).Please note, additionally, that j Γ = 0 at Γ W due to the no-slip condition.

Objective 2: Scalar shear stress metric
The target here is to minimize a scalar shear stress metric which refers to the squared volume-specific shear stress defined in an observation region of the fluid domain (Ω o ).The computation of the objective functional follows from Similar objective functionals have previously been considered for the shape 23 or topology optimization 34 of arterial bypass grafts.However, previous works mostly focused on minimizing the squared norm of the strain-rate tensor (S ik S ik ), as it relates to mechanical blood damage (hemolysis).Here, we augment the squared norm by the effective viscosity of the fluid.The monitored quantity thereby relates directly to shear stresses, which is the cause of hemolysis.Mind that viscosity acts as a constant scaling factor on the stress for a Newtonian fluid, while it augments the non-linear dependency of the stress on the strain-rate for a non-Newtonian fluid.Therefore, a stress-based objective facilitates a direct comparison between Newtonian and non-Newtonian adjoint methods.Due to the clinical relevance of this objective and practical restrictions that might arise in this context, for the minimization of this objective we target shape alterations that don't coincide with the boundaries of the observation region.The latter are also considered not to coincide with the domain's inlet (Γ in ) and outlet (Γ out ). Formally, we thus consider where Γ o denotes the boundaries of the observation region.It then holds, Even though the objective is defined in the domain, boundary contributions arise as shown in Equations ( 24).We thus have contributions that enter the adjoint system in sections of both the domain and its boundary.A detailed analysis for the derivation of Equations ( 24) as well as their impact on the adjoint system can be found in Appendix B.

Boundary conditions
The boundary conditions (BCs) for closing the generalized adjoint equation system are deduced by demanding that all surface integrals in Equations ( 14)-( 16) vanish.The aforementioned integrals may vanish either by demanding that the bracketed quantities, involving adjoint quantities, are set to zero or naturally by the BCs of the primal problem.

Inlet and wall
Primal velocity is usually fixed on Γ in and Γ W and we can thus claim u i = 0.The adjoint velocity can therefore be deduced in these patches by demanding on Γ in and Γ W . Based on the above, we can prescribe BCs for the adjoint velocity components on the inlet and wall patches.
As for the adjoint pressure, a lack of a constraint to be satisfied leads us to impose a zero Neumann boundary condition to assure the well posedness of the adjoint equation system, see Reference 26.

Outlet
A Dirichlet as well as a zero Neumann boundary condition are prescribed for pressure and velocity, respectively, on Γ out .Therefore, we claim leading to the identical cancellation of the surface integral in Equation ( 14) and the second surfance integral in Equation (15).The constraints that we need to impose on Γ out , thus, reduce to TA B L E 3 Boundary conditions for closing the primal and adjoint equation systems.
û n n = 0 , û t i = 0 Equation ( 33) For both investigated objective functional, j Γ  μ = 0 and Equation ( 28) can be developed as by decomposing the remaining term into normal and tangential components.In specific, we denote the normal component of a vector by subscripting it with n while the tangential ones with the superscript t.A standalone n refers to the spatial normal direction.We assume that, in the converged state of the cases studied in this work, u n n vanishes on the outlet and thus we define û t i = 0 to satisfy Equation ( 29).Subsequently, we decompose Equation ( 27) into its normal and tangential components, namely, By further developing term T1 of Equation (30), it is shown that since we already imposed û t i = 0 on the outlet.Furthermore, term T2 also vanishes owing to X ij ∼ S ij .Equation ( 30) thus reduces to From Equation ( 33) one can prescribe a Dirichlet BC for adjoint pressure.As regards the normal adjoint velocity component, we impose a BC based on the satisfaction of adjoint continuity (cf.Equation 17).We acknowledge that by doing so, Equation ( 31) is not necessarily satisfied, however, our computations showed that the satisfaction of the latter does not significantly impact the computation of the shape sensitivity.In the above, we intentionally omitted reference to BCs of primal and adjoint viscosity.This is because their computation is based on algebraic relations and thus no explicit boundary conditions are required.However, boundary values are required for the computation of related expressions, for example, Equation (33).Therefore, in all investigated cases, for the value of primal and adjoint viscosity we employ a zero gradient extrapolation to boundary patches from the nearest cell, which can be loosely translated to a Neumann boundary condition.The employed primal and adjoint BCs for all investigated cases are collected and presented in Table 3.

Shape sensitivity
The directional derivative of the Lagrangian in the direction of the control reads Furthermore, utilizing the fact that the total variation of the residual governing equations also vanishes, it holds, By substitution of Equation ( 36) in (34), It is easy to notice that the integral of Equation ( 37) corresponds to ∑ 5 k=1 ( y I k ⋅ y) of Equation ( 9), that has already been computed for the derivation of the adjoint equations and BCs.In specific, this term amounts to the integrals derived in Equations ( 14)-( 16), excluding objective contributions.Furthermore, the objective functional considered in this work are defined in sections of the domain that are not free for design and we thus consider  c J ⋅ c = 0. Based on the above, we can expand Equation ( 37) and write, As regards the volume integrals of Equation (38), one notices that they vanish by virtue of the adjoint equations, in case of zero or negligible * domain objective contributions, as is the case for the investigated objective functional of this article.Furthermore, we note that  c  ⋅ c only relates to Γ D ⊂ Γ W and since we have already imposed û i = 0 on Γ W (see Table 3), Equation (38) simplifies to Here, we follow the methodology introduced in Reference 36 and frequently employed thereafter 21,26,28 and approximate u i on the design wall through a first order Taylor series expansion of the no-slip condition, namely, where (c k n k ) illustrates that any changes to the shape are only plausible in the normal direction of the design surface.
Changing to index notation on the LHS of Equation ( 39) and utilizing Equation (40) results in We assume that continuity holds also along the boundaries and neglect curvature terms, thus u n n = 0.This leads to Based on the notation in Reference 37 we define the scalar shape sensitivity distribution for both objective functional as which we can then compute in each discrete grid face of our design section and guide our shapes toward improved states.
As seen in ( 43), the frequently employed sensitivity expression of a Newtonian fluid is augmented by a non-Newtonian contribution.

Numerical method
The numerical solution of the primal and adjoint equations is based upon the finite-volume method (FVM) employed by FreSCo + . 38The segregated algorithm uses a cell-centered, collocated storage arrangement for all transport properties.Possible parallelization is realized by means of a domain decomposition approach. 39,40Both primal and adjoint pressure-velocity coupling is based on the SIMPLE method.The convective fluxes in primal equations are discretized through the QUICK scheme.Due to the negative sign of the convection term in the adjoint momentum, a downwind analogy to QUICK is used.A detailed explanation of the discretized numerical solution approach for both primal and adjoint equations can be found in Reference 41.The primal [adjoint] non-Newtonian viscosity, μ [ μ], is computed at every cell center, on each SIMPLE iteration, based on the primal [adjoint] velocity of the previous iteration.The values of primal and adjoint viscosities are then obtained on the boundaries of the domain by means of a zero order extrapolation.
As regards the adjoint momentum equations for objective 2 (see Appendix B), additional contributions stemming from the objective are added explicitly to the RHS of the linearized equations, in sections of the computational grid embedded in the observation region.The shape sensitivity distribution (Equation 43) is evaluated after convergence of adjoint equations as a post-processing step.

VALIDATION STUDIES
This section showcases the validity of the previously derived surface sensitivity outlined by Equation (43).To this end, the computed shape derivative is compared against locally evaluated second order FD in a 2D exemplary geometry, see Figure 1A, and we consider solely objective 1.For the purposes of this section, similar geometries could be used, for example, relating to axisymmetric flows, however, they are skipped from this article for the sake of brevity.The geometry under investigation relates to a generic arterial stenosis, 42 the shape of which is modeled through a cosine curve, namely, where x 1 = L∕2 denotes the longitudinal position of the stenosis center, x 0 refers to half the length of the stenosis and  = H∕4 to its maximum height.The geometry is discretized with 12000 control volumes (CVs).The employed structured grid is progressively refined toward the design section, as shown in Figure 1B.In all investigated cases, boundary conditions follow from Table 3, where u in i corresponds to a parabolic, Poiseuille velocity profile, namely, We consider two different mean velocities (U) resulting to Re = UH∕ μN = 1.2 and Re = 12, where μN corresponds to the Newtonian viscosity.Our studies are confined to these low Reynolds numbers to ensure that the resulting shear rates are in regions of non-Newtonian interest for the investigated cases in this section.In the circulatory system the Reynolds number can vary from 0.001 (in a capillary) to 7500 (peak value for the aorta). 43n a simple, almost uni-directional shear flow, as the one considered herein, the Reynolds number is well correlated to the range of resulting shear rates and thus viscosities.However, in more involved cases, as the bypass graft studied in the application Section 4, the range of flow patterns increases and the flow might exhibit larger recirculation regions, localized strong accelerations and so forth.This also expands the range of shear rates, and therefore non-Newtonian properties might be relevant for the hemodynamics even in cases with high inlet Reynolds number.This applies in particular to situations where the natural blood flow is disrupted, either by necessary interventions or pathological factors, giving rise to unnatural flow patterns.
The FD study is realized by perturbing a finite amount of discrete faces on Γ D in their normal direction by ±.The objective functional is then computed for each primal simulation of perturbed shape and the local surface sensitivity is estimated for a surface patch k through where c k i represents the position of the design face centre,  is the magnitude of the perturbation, n k i is the normal vector at c k i and A k is the corresponding perturbed area.The latter is used as a support area for unit consistency w.r.t.Equation (43).
Overall, for the discrete estimation of surface sensitivity through the FD method, one needs to perform 2N k -where N k is the total number of perturbed faces on Γ D -primal simulations while on the other hand only one primal and one adjoint simulation is required for the continuous approximation of sensitivity using the adjoint method.This is one of the main advantages of the latter.For a reliable comparison of FD-based sensitivities against adjoint-based ones, we first study the impact of the perturbation magnitude on the former.Figure 2 shows that the employed sizes lie satisfactorily within a regime dominated by a linear response for both considered Re numbers.Therefore, for all subsequent FD-based results a perturbation magnitude of 10 −5  is employed.Figures 3 and 4 show FD-based (46) and adjoint-based (43) shape sensitivity estimations for Re = 1.2 and Re = 12, respectively.In specific, computations are realized for both Newtonian and non-Newtonian fluid models.For the latter, we further consider a reduced non-Newtonian model in which the primal problem is modeled through the respective non-Newtonian description while the dual follows from the Newtonian set of adjoint equations using the primal viscosities, that is, derivatives w.r.t. the "frozen" viscosity are neglected.As shown, the FD-based results (S k FD ) satisfactorily match their adjoint-based counterparts (S adj ) for both Newtonian and non-Newtonian fluids and for both considered Re numbers.However, when the reduced non-Newtonian model is considered, the adjoint-based shape sensitivity estimation ( Ŝadj ) deviates from the respective FD results.To quantify the significance of this deviation, for each considered Re and non-Newtonian model, we calculate T A B L E 4 Estimation of the relative difference between the shape sensitivities predicted from FD and the reduced ("frozen") non-Newtonian approach Ŝadj according to Equation (47).

Case
where N m is the number of perturbed faces in the most sensitive region (x − x 1 )∕x ∈ [−0.034, 0.034] and Ŝi adj denotes the local, reduced non-Newtonian, adjoint-based shape sensitivity at the same face.Results are collected and presented in Table 4.For Re = 1.2, the largest deviation is found for the Carreau model followed by the Power-Law.For Re = 12 the error of the former is significantly decreased while the latter's remains relatively unchanged.Furthermore, we note that the modified-Casson model produces the lowest error for both considered cases and that an increase in the Re number closes the gap between the estimated values for all non-Newtonian models.
The reported results are in accordance to our expectations.As the Reynolds number increases, the non-Newtonian effects fade out and the fluid can satisfactorily be described by a constant Newtonian viscosity.Therefore, the error between the FD-based results and those of the reduced non-Newtonian adjoint model, which essentially corresponds to the classic Newtonian adjoint model, decreases.Figure 5 shows the effective viscosity of each model w.r.t. the shear rate.As can be seen, the primal viscosity lies in a close neighborhood of the corresponding Newtonian one, for the modified Casson model, while it deviates quite significantly from the Newtonian value for Power-law and Carreau models in the low shear-rate regime.The aforementioned observation explains the increased errors of the (PL) and (C) models in comparison to the (MC) model displayed in Figures 3 and 4. As regards the extent to which the error drops by increasing the Reynolds number, it is attributed to the nature of each model.In the cases studied herein, the maximum computed shear rates were approximately 50 and 500 (1/s) for Re = 1.2 and Re = 12, respectively.In this range, the (C) and (MC) models start asymptotically converging toward the (N) viscosity thus justifying their small error for Re = 12.On the other hand, in the investigated range of shear rates, the (PL) model results in effective viscosities that equidistantly differ from the Newtonian one.We attribute the almost constant sensitivity difference displayed by the (PL) in Table 4 to this observation.It is noted that our observations consider a direct relation between viscous diffusion and shape sensitivity, which does not naturally hold.That is, we assume that the difference between the FD-based sensitivity results and the reduced adjoint results exclusively depends on the difference between Newtonian and non-Newtonian effective viscosity.While there is indeed a direct dependence between the adjoint-based sensitivity outlined in Equation ( 43) and the effective viscosity, the expression is far more involved than that.However, due to the low Reynolds numbers considered, momentum transported is governed by viscous diffusion.This primal flow behavior naturally reflects on the adjoint problem and thus the shape sensitivity estimation closely relates to the nature of the effective viscosity of the primal problem.
Finally, we would like to remark on the trend of S adj and Ŝadj .While the differences can be quite significant, especially in regions where non-Newtonian properties are dominant, the overall trend of the two curves is fairly similar, that is, one could approximate Ŝadj ≈ kS adj , where k represents a constant value.Naturally, we cannot have an a priori knowledge of what k should be, however, in the context of a gradient-descent optimization approach, the difference between a reduced and a complete non-Newtonian model would essentially boil down to a step-size issue.Nevertheless, while this observation holds for the cases studied in this section, it does not necessarily hold for more involved cases.

EXEMPLARY APPLICATION
In this section we employ a shape optimization procedure, on a 3D idealized arterial bypass-graft, for the minimization of the scalar shear stress metric, defined as objective 2 (see Section 2.1).The initial geometry of the investigated case is shown in Figure 6.The study of this configuration is motivated by its clinical relevance in surgical procedures.Bypass-grafts are often used to divert blood around narrowed or occluded parts of an artery.For the cases studied herein, the bottom leftmost part of the impaired artery is modeled as occluded.The success of the surgery has been shown to strongly depend on the hemodynamics, 44 partially resulting from the shape of the implanted graft.The use of objective 2 as our cost functional is motivated by the significance of shear stresses in the flow field when designing shapes associated with blood flow.Excessively high stresses induced by peculiarities of the blood flow may damage the red blood cells or even lead to thrombus formation. 21,23Furthermore, the objective functional is constructed in such a way that the observation region is detached from the design.This allows us to study the impact of upstream shape changes of the graft-where the surgeon can act before or during the operation-on downstream hemodynamics.In all subsequent studies, we consider a complete non-Newtonian (CNN) adjoint model and a reduced (RNN) one for the estimation of the shape sensitivity.On the latter, derivatives of viscosity are omitted and μ is considered "frozen" for the dual problem.

Investigated cases
In all investigated scenarios, blood's density, , is set to 1056 kg∕m 3 while μ follows from the Carreau model.Figure 7A shows the side view of the initial geometry, annotated with key geometrical parameters.In specific, we consider x = 0 to coincide with the inlet patch and the total longitudinal length L, from inlet to outlet, to be 93 mm.   3.

Reynolds number (Re) Computational grid (i) Number of CVs
To ensure that the evaluation of our objective functional J is independent of the spatial discretization, a grid study is conducted, as presented in Table 5.The numerical evaluation of J, see Equation (23), follows from the mid-point approximation rule.The employed structured numerical grid is displayed in Figure 7B and consists of approximately 300k CVs.It corresponds to the first level of refinement (second grid), since no significant change in the estimation of the objective is reported for the next level, as shown in Table 5.

Shape optimization procedure
The shape optimization studies are realized in terms of a CAD-free optimisation approach, in which the computational grid is successively adjusted after each primal and adjoint run.Once the shape sensitivity is computed, an auxiliary problem is solved to determine the shape gradient from which the discrete displacements of boundary and internal nodes can be estimated based on a gradient-descent method.In this article, we employ a Steklov-Poincaré type of approach to approximate the descent direction (gradient) g i , which reads where q = 1∕(W D + )-with W D denoting the wall normal distance and  a small value to avoid numerical error on the boundary-is used to propagate g i smoothly in the domain.This approach has been successfully used in recent large-scale engineering applications, see for example, References 48 and 49.We also point the interested reader to Reference 50 for a comprehensive comparison, from a mathematical basis, of different approaches.Equations ( 48)-( 50) are solved in a FV fashion, from which g i is computed for each cell center.Subsequently, the values are mapped to the nodal positions x n i using an inverse distance weighting, also known as Shepard's interpolation. 51Here, we use superscript n and c to denote nodal values and cell center values, respectively where N c denotes the total number of adjacent cells to node n.The grid is then updated by applying where  denotes a constant step size, and geometric quantities are recomputed for each CV.The topological relationships of the grid remain unchanged and the optimization process continues by restarting from the previous optimization step.Due to the employed iterative optimization algorithm and comparably small step sizes, field solutions of two consecutive shapes are usually nearby and thus the computational time of a complete optimization process is significantly reduced compared to starting each optimization step from scratch.The above process is repeated until convergence or until the maximum number of shapes is reached.The optimization process is outlined in Algorithm 1.In all our investigations, the step size is calculated in the first optimization iteration as and remains constant throughout the optimization process.This allows us to further control the shape updates and avoid large deformations that might break the computational grid, by selecting the maximum cell displacement of the first optimization iteration through g max , which is set to g max = D in ∕200, for all cases studied herein.The step-size identification (53) helps to control step-size issues when comparing the RNN and CNN approaches, compare, Section 3. models are mainly located close to the end of the design section and notably differ between the two approaches.The difference is even further highlighted when the sensitivity is passed on the auxiliary problem to compute the deformation in Equations ( 48)- (50).The resulting displacement field is significantly different, especially in the cuff area of the graft.We note, therefore, that neglecting contributions from viscosity derivatives on the adjoint side results in differences on the resulting sensitivity maps and the corresponding displacement fields.However, when the input velocity is increased we note that these differences are visibly decreased as shown in Figures 10 and 11.In specific, we note that while for Re in 3 there are still some notable differences in the sides of the cuff, these differences almost completely vanish for the highest considered Reynolds number.These observations are in accordance with our validation studies in Section 3. We would also like to point the attention of the reader to the maximum value of displacement applied on the first optimization iteration, which is the same between all three considered cases.This is due to the applied step size scheme as outlined in Equation (53).Based on the the observations made in Section 3, one would expect that the resulting displacement fields would be, if not the same, very similar.We show here that this is not the case for lower Reynolds numbers where notable differences occur.

Results and discussion
In order to further quantify the impact of considering the CNN adjoint model, we perform three (one for each considered Reynolds number) additional identical primal simulations on the initial shape, using the Newtonian model.Figure 12 shows the percentage difference between the computed objective value using the non-Newtonian Carreau model and the Newtonian model.Figure 13 shows the viscosity field for each considered Re in .As shown, the higher the inflow flux the less prominent the primal non-Newtonian properties.The leftmost part of high effective viscosity values corresponds to the occluded section and is of little importance.Compared to the Newtonian viscosity, elevated values are observed near separation regions in all cases.However, as the influx is increased, the flow is stabilized and separation or recirculation regions fade out.The extent to which the objective computation is influenced by the use of a non-Newtonian model is also reflected by the impact of using the CNN adjoint model for the computation of the shape sensitivity.Furthermore, it is interesting to note how these differences may affect the final result of the optimization.Figure 14 shows the relative decrease of the objective compared to the initial value throughout the optimization process.We note that optimizations for Re in 1 met the convergence criterion, while those for Re in 2 and Re in 3 reached the maximum number of optimization steps.Overall, the presented results meet our expectations and conform with those presented in Section 3.For the smallest Re number the difference between the relative decrease of the objective function for the CNN and RNN models is approximately 3%, while for the higher Re in 1 and Re in 2 the difference is negligible.This shows us that for low enough Re numbers, where the shear-rates are sufficiently small to activate the shear-thinning properties of blood, the CNN model is not only more accurate on the estimation of sensitivity (see Section 3) but also more efficient in finding a sufficiently optimized shape compared to RNN.For larger Re numbers, however, we note an equivalence of the two methods.Furthermore, we note the difference in percentage decrease between the three cases.This is attributed to the nature of the objective functional which amounts to subsequently smaller values for flows of low shear-rates and thus the influence of the shape on it becomes of less importance.
In Figure 15 we show the differences of the optimized shapes produced by the CNN and RNN models.We note that the difference between the two optimized shapes reduces as the Reynolds number increases.However, interestingly enough, the distribution of shape differences changes significantly for each individual considered case.This shows us that even though the objective reached similar levels of reduction between the two adjoint models, different shapes can be expected.This is justifiable due to localized areas of non-Newtonian interest that might occur even in higher Reynolds number cases, for example, recirculation regions.Therefore, while for the cases studied herein, the complete optimization history seems to be mostly driven by the Newtonian nature of the adjoint equations, if one neglects non-Newtonian contributions, some details might be overlooked upon.
Furthermore, it is interesting to note the sudden cross-sectional change from the design to the non-design section for the optimized shapes obtained for Re in 1 .This is due to the construction of the shape deformation problem, which follows from Equations ( 48)- (50), as well as the resulting shape sensitivity, as shown in Figure 9, which is initially condensed in the region close to the non-design section.This issue can be circumvented by applying an additional filter to the descent direction field g i , that is, in the neighborhood of the transition between Γ D and Γ ⧵ Γ D , so that shape and mesh deformations smoothly vanish close to a non-design section, compare, References 52 or 53, where such techniques are applied for a 2D and a 3D case.At the interest of outlining the differences between the RNN and CNN adjoint models, these techniques are not applied herein, as they would alter the solution of the optimization problem.
Additionally, we would like to discuss the optimized shape obtained for Re in 2 .One would expect that the kink created in the cuff region would deteriorate the hemodynamics and lead to higher stresses, which is utterly what is observed by objective 2. This is indeed the case upstream of the flow, where the quantity monitored by the objective has increased by 26% from initial to optimized case, due to a recirculation zone appearing downstream near the kink.However, in Figure 16 we show that the optimizer led to an artificial "nozzle-like" shape that smoothens the flow in the observation region leading to a decrease of the objective functional, which is only measured in the observation region, by 16%.It is important, therefore, that in the context of biomedical applications, the decision of the objective is informed by clinical knowledge.
Table 6 shows the wall clock time required for each optimization process.The notable increase in wall clock time for the cases corresponding to Re in 1 is due to the additional computational steps required for convergence of the CNN approach.For Re in 2 and Re in 3 the maximum number of iterations were reached and we note a small difference between the two models.This is attributed to the additional computations necessary for the CNN model.Note: Each primal/adjoint simulation was parallelized using 96 cascade 9242 CPUs.Final column denotes the percentage increase of wall clock time when using the CNN compared to RNN models.
Overall, the optimization studies conducted in this section confirm our hypothesis as regards the impact of neglecting viscosity derivatives in the adjoint problem.While the impact may be subtle for flow scenarios of high shear rates, it becomes noteworthy for low shear rates.

Unsteady primal simulations
Due to conducting the optimization studies in steady-state for three fixed inlet conditions, we also arrive at three different optimal shape variants.In order to optimize the initial shape for the complete heartbeat, one must consider an unsteady adjoint optimization run.While such unsteady adjoint approaches go beyond our scope, we evaluate the optimized shapes in an unsteady framework as an a posteriori step.In specific, the initial as well as the optimized shapes, produced from the CNN adjoint method, are simulated for two heartbeats.A Womersley velocity profile 54 based on the volumetric blood flux, shown in Figure 8, is prescribed at the inlet, while a zero pressure condition was employed along the outlet.The blood properties agree with the properties used in the optimization runs.
Figure 17 shows a comparison of the temporal evolution of the objective function for the second objective.The figure depicts the performance of the three optimized shapes (black line), that is, (A) Re in 1 , (B) Re in 2 , and (C) Re in 3 , in comparison to the initial configuration (red line).Figure 18 displays the temporal evolution of the corresponding relative objective function change, which portraits instantaneous improvements (negative J opt − J ini ) and deteriorations (positive J opt − J ini ) in comparison with the initial shape.The shape produced by optimization at Re in 1 shows small absolute improvements for the complete heart-cycle.However, these improvements amplify in a relative sense for time instants when the influx is low and thus closer to the conditions on which the steady-state optimization has been performed.The shapes produced by Re in 2 and Re in 3 show a more significant absolute improvement in periods of increased influx.When the influx is lower, the performance of the optimized shape can deteriorate as compared to the initial shape.
This study shows, that the time-averaged transient performance of all shapes optimized under steady conditions exhibits a beneficial response in comparison to the initial shape.However, the improvement is temporally localized around a time integral of the steady optimization conditions and can even reverse for time instances of vastly different influxes.Therefore and in order to perform a clinically relevant design, the need for an unsteady adjoint-based shape optimization is highlighted.

CONCLUSIONS AND OUTLOOK
We have proposed a system of adjoint equations, dual to the Navier-Stokes equations with a generalized constitutive model for non-Newtonian fluids.The adjoint model is derived for three frequently used non-Newtonian models, in the context of blood flows, and presented in a generalized form.We further specified the model for two objective functionals, one commonly studied in CFD-based shape optimization and one relating to blood-damaging factors with clinical relevance.Adjoint boundary conditions, closing the dual problem, and an expression to approximate the shape sensitivity are proposed.The numerical solution of both primal and adjoint problems is based on the finite-volume method.The proposed model is validated against second-order accurate finite-difference results on a 2D exemplary geometry loosely related to an arterial stenosis.A comparison is made with respect to a reduced adjoint model, where the primal is modeled through the use of a generalized constitutive model while derivatives of viscosity are omitted in the adjoint equations.Finally, we included the model in a CAD-free, gradient-based shape optimization framework and applied it for the minimization of a shear stress metric on an idealized 3D arterial bypass-graft.The considered metric is used to indirectly quantify the potential of the implanted graft to damage the diverted blood.The goal of this study was to minimize this potential by altering the shape of the graft.We have considered three physiologically relevant blood fluxes for a steady-state optimization of the shape from which we identified different optimized shapes.An a posteriori evaluation of the shapes in unsteady conditions showed that the beneficial behavior of each shape is localized around the time for which the design was performed.
It was shown that when an objective functional, defined exclusively on the boundaries of the domain, is selected, the sole difference between the traditional adjoint Navier-Stokes equations and the ones proposed herein refers to an additional contribution in the adjoint momentum equations.The proposed sensitivity derivative was also enhanced with an additional term relating to the adjoint viscosity.Through the validation studies, we have observed that neglecting these contributions may result to inaccurate shape sensitivities.Further, we showed that its impact increases as the Reynolds number decreases, a behavior that we attribute to the shear-thinning properties of the investigated models.This was also verified by our findings on the shape optimization studies.In specific, we showed that considering a reduced non-Newtonian adjoint model may result to differences of up to 3% on objective reduction as well as different selected shapes.In terms of numerics, the implementation of the proposed model is analogous to its primal counterpart and the necessary adjoint viscosity is computed after each adjoint velocity-pressure coupling iteration.We showed that in the context of an optimization process the additional computational effort does not lead to more than 10% increase in wall clock time for the same optimization iterations.
As a bottleneck of our optimization studies, we identified the competitive behavior between the scalar shear stress metric and the shear-thinning behavior of blood.For flows of low shear rate, while the non-Newtonian properties are activated, the objective becomes less significant by construction and vice versa.To this extent, we are also interested in applying the proposed methodology for shear-thickening fluids, in which we expect a reversal of the trends presented in this article.Furthermore, we recognize that by considering steady-state flow simulations, the optimized shapes do not necessarily represent beneficial configurations for the complete pulsatile blood flow, as also shown from our a posteriori evaluation studies.
We are working toward extending the model to involve unsteady simulations as well as to incorporate it in a fluid-structure interaction framework.This would enable us to simulate more realistic flow scenarios of clinical relevance.Both methods require efficient strategies for handling the extensive amount of primal data due to the reversed time propagation of the adjoint.We point the reader to Reference 55 for a recently proposed such strategy.Finally, we target to extend to a robust shape optimization framework.In this context, the fluid's constitutive model parameters can be applied with a certain degree of uncertainty and the impact of this uncertainty propagated to the final result.
Here S ij = 0.5(u i ∕x j + u j ∕x i ) and higher order terms (HoT) are neglected.We define X ij = 2k(n − 1) ̇(n−3) S ij , on which it holds X ij ∼ S ij − −−−− → X ij = X ji .Equation (A2) is further expanded as

A.3 Modified-Casson model (MC)
In the case of a Modified-Casson fluid Equation (A1) can be written as S ij and the same methodology as described in Section A.2 is followed, then Equation ( 12) is reached.

A.4 Carreau model (C)
In the case of a Carreau fluid Equation (A1) can be written as ) n−3 2 Λ 2 S ij and the same methodology as described in Section A.2 is followed, then Equation ( 12) is reached.

APPENDIX B. ANALYSIS RELATED TO OBJECTIVE 2
For the analysis herein we consider an illustrative 2D geometry as shown in Figure B1.We consider an observation region Ω o ⊂ Ω with boundaries Ω o = Γ o .As shown in Figure B1, (Γ D ∪ Γ in ∪ Γ out ) ∩ Γ o = ∅ applies.
By perturbing the objective,

F
I G U R E 1 (A) Sketch of the 2D geometry for the FD study.Dashed line corresponds to the design section (Γ D ).The origin coincides with the leftmost lowest point of the geometry.(B) Detail of the employed numerical grid near the design section.

2 3 F I G U R E 4
Influence of perturbation magnitude  on S k FD for c k i = (x 1 , ).The computed objective functional of the unperturbed and perturbed shapes are denoted by J 0 and J * , respectively.Investigations for Newtonian fluid at (A) Re = 1.2 and (B) Re = 12.[Colour figure can be viewed at wileyonlinelibrary.com]Study for Re = 1.2:Shape sensitivity estimations for a (A) Power-Law, (B) Modified-Casson, and (C) Carreau fluid.Black (continuous), red (dashed), and blue (densely dotted) lines correspond to Newtonian, non-Newtonian, and reduced non-Newtonian adjoint-based estimations, respectively.Circles and triangles denote discrete FD results for Newtonian and non-Newtonian fluid, accordingly.Perturbation magnitude ∕ = 10 −5 .[Colour figure can be viewed at wileyonlinelibrary.com]Study for Re = 12: Caption and legends as in Figure 3. [Colour figure can be viewed at wileyonlinelibrary.com]

5
Effective viscosity for (A) Power-Law, (B) Modified-Casson, and (C) Carreau viscosity models.The Newtonian's fluid viscosity is included in black [continuous line] for comparison purposes.Axes in logarithmic scale.[Colour figure can be viewed at wileyonlinelibrary.com]

8
For the inlet and outlet diameters, it holds that D in = L∕6.5 and D out = L∕10, respectively.The offset of the proximal end from the inlet is  = L∕11.6 and the size of the cuff, that is, connection between the graft and the impaired artery, is l = L∕3.3.In our optimization studies we consider Γ D = {(x, y, z) T ∈ Γ W ∶ x ∈ [L∕6, 2L∕3]}, thus including in our design space the graft and the proximal end, areas in which the surgeon can act upon.The observation region is defined downstream asΩ o = {(x, y, z) T ∈ Ω ∶ x ∈ [L∕1.4,L∕1.05]}.As concerns the simulation of similar problems, one would usually consider a transient flow or even fluid-structure interaction, to accurately resolve the flow phenomena, see for example, References 45 and 46.While this goes beyond the scope of this article, we are still interested in optimizing the given geometry under physiological conditions and thus we employ the volumetric flux, shown in Figure8, as a reference point for the prescription of BCs.The flow pulse was extracted from Reference 47 and scaled down so that physiological velocities arise for the given geometry.As shown in Figure8, we select three distinct time instants for our steady-state adjoint-based shape optimization Γ Γ Γ Γ F I G U R E 6 Initial geometry of the idealized arterial bypass-graft.U R E 7 (A) Side view of the initial geometry.(B) Perspective view of the utilized computational grid.Physiological volumetric blood flux over a heartbeat (≈ 0.8s), adjusted from Reference 47. Circles denote the selected fluxes for the optimization studies herein.[Colour figure can be viewed at wileyonlinelibrary.com]TA B L E 5 Results of the mesh dependence study.Index i refers to the mesh level.Refinement is realized in all directions of the investigated O-H-grid.

Figures 9 -F I G U R E 9
Figures 9-11 show the sensitivity and the resulting displacement maps of the design section (Γ D ) for the first optimization step.As can be seen for Re in 1 in Figure9A,B, the computed sensitivities resulting from the RNN and the CNN adjoint

F 11 12
I G U R E 10 Study for the highest mass flux at Re in Study for the medium mass flux at Re in 3 = 295: Caption as in Figure 9. [Colour figure can be viewed at wileyonlinelibrary.com]Objective percentage difference between the objective computed on the initial shape with a Newtonian and a non-Newtonian model.Computation follows from |V1 − V2| ⋅ 2 ⋅ 100∕(V1 + V2), where V1, V2 denote the objective value as computed using the Newtonian and non-Newtonian model, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

13 14
Effective viscosity field for (A) Re in 1 , (B) Re in 2 , and (C) Re in 3 .Computations on the initial shape.[Colour figure can be viewed at wileyonlinelibrary.com]Relative decrease of objective J during the optimization run for Re in = (A) 63, (B) 507, and (C) 295, using the complete (cross marks) and reduced (circle marks) adjoint models.The objective decrease is computed as (J i − J 0 )∕J 0 ⋅ 100%, where the index denotes the number of shape.[Colour figure can be viewed at wileyonlinelibrary.com]

15
Differences between the optimized shapes using the CNN and RNN adjoint models on the z − x plane.From top to bottom results correspond to Re in 1 , Re in 2 and Re in 3 , respectively.Left: Outlines of the initial (blue line) and the optimized design section by employing the CNN (red line) and RNN (black line) models.Right: Normalized cell coordinate distance contours plotted on the optimized shape produced by the RNN model.The computation of this field follows from ||x CNN − x RNN ||∕D in , where x CNN and x RNN denote the cell coordinates of the optimized shapes produced by the CNN and RNN models, respectively and D in corresponds to the diameter of the inlet.[Colour figure can be viewed at wileyonlinelibrary.com]

F I G U R E 16
Study for Re in 2 : Comparison of velocity profiles in the observation region for initial (left) and optimized (right) shape.Cross sections correspond to the start, middle and end of the observation region.Figure augmented by a detail of the designed region in the vicinity of the observation region.Optimal shape corresponds to the one obtained from the RNN approach due to similarity with the one obtained with the CNN approach.[Colour figure can be viewed at wileyonlinelibrary.com]

17 18
Evolution of normalized objective J∕J ref for initial and optimized shapes in steady-state for (A) Re in 1 , (B) Re in 2 , and (C) Re in 3 , during two heartbeats.Normalization based on J ref = 100 ( μN u max ∕D in ) 2 V o , where μN , u max and V o are the Newtonian viscosity, maximum velocity magnitude at the inlet during a heartbeat, and the volume of the observation region, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]Evolution of the relative change between the objective computed at the initial (J ini ) and at the optimized shapes (J opt ) for (A) Re in 1 , (B) Re in 2 , and (C) Re in 3 , during two heartbeats.[Colour figure can be viewed at wileyonlinelibrary.com] U R E B1 Exemplary domain for the analysis of objective 2. Observation region (Ω o ) is shown with dashed line while the design regions (Γ D ) are shown in dashed-dot lines.
2 to cover a broad range of Reynolds numbers.The latter is calculated as Re in = UD in ∕ μN , where μN corresponds to the Newtonian viscosity and U denotes the mean of the parabolic inlet velocity, with Q=U D2 studies

TA B L E 6
Measured wall clock times (t wc tot ) for a complete optimization process for each considered case.