Weakly imposed boundary conditions for shear-rate dependent non-Newtonian fluids: application to cardiovascular flows

This paper presents a stabilized formulation for the generalized Navier-Stokes equations for weak enforcement of essential boundary conditions. The non-Newtonian behavior of blood is modeled via shear-rate dependent constitutive equations. The boundary terms for weak enforcement of Dirichlet boundary conditions are derived via locally resolving the fine-scale variational equation facilitated by the Variational Multiscale (VMS) framework. The proposed method reproduces the consistency and stabilization terms that are present in the Nitsche type approaches. In addition, for the shear-rate fluids, two more boundary terms appear. One of these terms is the viscosity-derivative term and is a function of the shear-rate, while the other term is a zeroth-order term. These terms play an important role in attaining optimal convergence rates for the velocity and pressure fields in the norms considered. A most significant contribution is the form of the stabilization tensors that are also variationally derived. Employing edge functions the edge stabilization tensor is numerically evaluated, and it adaptively adjusts itself to the magnitude of the boundary residual. The resulting formulation is variationally consistent and the weakly imposed no-slip boundary condition leads to higher accuracy of the spatial gradients for coarse boundary-layer meshes when compared with the traditional strongly imposed boundary conditions. This feature of the present approach will be of significance in imposing interfacial continuity conditions across non-matching discretizations in blood-artery interaction problems. A set of test cases is presented to investigate the mathematical attributes of the method and a patient-specific case is presented to show its clinical relevance.


Introduction
Cardiovascular blood flow simulations of clinical relevance require non-Newtonian blood flow models and invariably involve complex patient-specific geometries with branching arterial trees. The quality of the computed solutions is affected by the precise description of the geometric models as well as the quality of the computational grids. The use of boundary layer meshes helps in accurately calculating flow quantities such as wall shear stress (WSS) and surface pressure at the arterial walls, and several mesh adaptation techniques have been proposed in the literature for local refinement along the boundaries [1][2][3]. This however comes at an increased computational cost due to the insertion of elements in the boundary layer region to achieve higher spatial resolution of the velocity and pressure gradients [1,3]. In this context a numerical method that can yield higher spatio-temporal accuracy of the gradients in the fields and can facilitate the coupling of the primary fields across the blood-tissue interaction surfaces via weakly enforced continuity equation for the velocity field can open the doors for developing clinically relevant computational techniques.
The methods for weakly enforced no-slip conditions embed the Dirichlet boundary conditions in the variational formulation rather than prescribing the value of the Dirichlet data directly at the nodal points. This strategy allows the fluid particles to slightly slide at the wall and it helps in capturing the high gradients of the velocity field near the boundaries. These methodologies are also reported to help increase the accuracy of the boundary layer fields [4]. Similarly, methods for weakly imposed boundary conditions have also been proposed for wall-bounded turbulent flows [5] where they behave like wall function models and facilitate better performance as compared to the strongly imposed boundary conditions. In addition, the notion of weak imposition of boundary conditions has been effectively employed in advection-dominated diffusion problems [6], in porous flows [7], and in buoyancy-driven flows [8]. These studies suggest the possibility of being able to use coarse discretizations in arterial flow analysis with dominant presence of the confining bounding surfaces, thereby significantly reducing the mesh sizes due to the computational efficiency engendered by the weakly imposed boundary conditions. Another advantage of the weakly imposed boundary conditions originates from the flexibility that these methods do not require nodally matched discretizations at the bloodartery interaction surfaces. The Dirichlet boundary is embedded in the variational equation via weakly imposed boundary conditions, thereby implicitly accounting for the continuity of the fields across the fluid-solid interface. Consequently, the computational meshes do not need to be nodally aligned, and this flexibility opens the door to develop methods for non-matching boundaries and interfaces [9][10][11] that are beneficial in complex geometries encountered in patient-specific models.
A literature review reveals that Nitsche-type methods have often been employed to variationally apply the no-slip boundary conditions [4,7,8,11,12]. The Nitsche method for constraining the primary unknown fields is comprised of the consistency terms and the penalty or stability term. A major technical issue has been the value of the stability parameter that largely remains unspecified by the theory. Since the condition number of the system is a function of the stability parameter, the optimal value of the parameter needs to be identified. This value needs to be small enough to preserve consistency of the method, while it needs to be sufficiently large for ensuring stability. If the parameter is too large, the stabilization term plays the role of penalizing the boundary conditions, and the formulation loses the variational consistency, thereby deteriorating the conditioning of the tangent matrix. If the value is too small, the formulation becomes unstable.
An important ingredient of the present paper is the use of shear-rate dependent non-Newtonian model for blood [13][14][15]. The flow of shear-thinning fluids generally gives rise to pseudo-plastic velocity profiles which is characterized by lower velocities due to a less convective flow but sharper and thinner boundary layers [15]. Therefore, the ability of the method to capture the boundary layer is of significant importance for the shear-thinning models than for the Newtonian models. The stabilized methods for this class of fluids were developed by the senior author, and the interested reader is referred to [15,[17][18][19][20][21].
This paper develops a stabilized formulation for the weak imposition of Dirichlet boundary conditions for non-Newtonian fluids. We derive the boundary terms by exploiting the edge-based functions that model fine scales in the Variational Multiscale (VMS) equations. Similar ideas have been applied to develop Discontinuous Galerkin (DG) methods for coupling multiple PDEs [22], in damage modeling for finite strain solid mechanics problems with weak and strong discontinuities across common interfaces [23], and in the immersed boundary conditions in fluid mechanics [10]. A unique contribution of the proposed method is the variationally derived consistency terms as well as the stabilization tensor that possesses a self-adjusting feature for optimal enforcement of the boundary conditions. This feature is highlighted with the help of numerical examples in Section 4.
An outline of the paper is as follows: Section 2 presents the governing equations and the constitutive model. Section 3 presents the derivation of the stabilized formulation. Section 4 presents numerical test cases to validate the proposed method and to investigate its mathematical and computational attributes. An application to patient-specific arterial geometry illustrates its clinical relevance. Conclusions are drawn in Section 5.

The strong form of governing equations
Let Ω ⊂ ℝ n sd be an open bounded region with piecewise smooth boundary Γ . The number of spatial dimensions, n sd , is equal to 2 or 3. The domain boundary assumes the usual split Γ = Γ g ∪ Γ ℎ and Γ g ∩ Γ ℎ = ∅, where Γ g and Γ ℎ are parts of the boundary with essential and natural boundary conditions, respectively. The governing equations and boundary conditions are given as: σ ⋅ n = h on Γ ℎ (4) where v and p are the fluid velocity and pressure, respectively. ρ is the density of the fluid, σ v is the deviatoric stress tensor, σ is the Cauchy stress tensor defined as σ = − pI + σ v , I is the identity tensor, f is the body force per unit mass, g is the prescribed velocity on the boundary Γ g , h is the prescribed traction on the boundary Γ ℎ , and n is the unit outward normal to the boundary of the domain Ω. Equations (1)-(4) represent balance of momentum, the continuity equation, the Dirichlet and Neumann boundary conditions, respectively. The shear-rate dependent deviatoric stress tensor is defined as where ε v is the rate-of-deformation tensor defined as ε v = 1 2 ∇v + ∇v T , γ is the shearrate defined as γ = 2ε v : ε v , and η γ is the effective viscosity which is a function of the shear-rate.

Non-Newtonian fluid models
In this work, we employ the power-law model and the Carreau-Yasuda model to represent the shear-thinning behavior of blood, which are the two commonly used models in computational hemodynamics [13][14][15]18,24]. The power-law model is: η γ˙= μγ˙n − 1 (6) where μ represents the viscosity of Newtonian fluids if n = 1. A drawback with this model is the singularity that appears at zero shear-rate, and various techniques have been proposed in the literature to establish a lower-bound on the effective viscosity.
In the Carreau-Yasuda model, the nonlinear function of the viscosity is given by η γ˙= μ ∞ + μ 0 − μ ∞ 1 + (λγ˙) a n − 1 /a where μ 0 and μ ∞ are asymptotic viscosities at zero and infinite shear-rate, respectively, and a, n and λ are empirically determined constitutive parameters. Coefficients a and n are non-dimensional parameters that control the shear-thinning or shear-thickening behavior of fluids in the non-Newtonian regime between the two asymptotic viscosities. The model reverts to the Newtonian fluid model by setting μ 0 = μ ∞ .

The mixed weak formulation
Our objective is to develop the formulation that weakly imposes the Dirichlet boundary conditions for the velocity field in Eq (3). We start from the classical mixed formulation and derive the expression for the Lagrange multiplier field at the Dirichlet boundary Γ g . Let w and q denote the weighting functions for the velocity and the pressure fields, and ψ denotes the weighting function for the boundary condition. The appropriate spaces of trial solutions and weighting functions are specified as follows: The mixed weak form associated with Eqs (1)-(4) is: Find v, p ∈ S × P, λ ∈ Q such that for all w, q ∈ V × P, ψ ∈ Q : (ψ, v − g) Γ g = 0 (13) where ⋅ , ⋅ Ω = ∫ Ω ⋅ dΩ is the L 2 Ω inner product. In this work shear-rate dependent deviatoric stress is considered. The mixed formulation in Eqs (12) and (13) consistently enforces the boundary condition Eq (3) through the Lagrange multiplier λ which acts as the numerical flux or the traction at Γ g . The downside of the mixed formulation is the presence of additional unknown fields and the instability associated with the inf-sup conditions. To overcome these issues while holding the virtue of variational consistency, we derive the expression for the Lagrange multiplier λ that will eliminate the auxiliary unknown field from the formulation.

The stabilized formulation for weakly imposed boundary conditions
This section presents the derivation of the stabilized formulation for the weakly imposed Dirichlet boundary conditions. We apply the VMS framework to the narrow band Ω ⊂ Ω along the Dirichlet boundary Γ g shown in Figure 1. Overall procedure underlying the derivation comprises three steps, (i) split the problem into coarse-and fine-scale subproblems, (ii) solve the fine-scale problem along the boundary Γ g , and (iii) embed the fine-scale models into the coarse-scale formulation.

Multiscale decomposition
We discretize the domain Ω into disjoint elements Ω e with element boundaries Γ e , such that Ω = ⋃ e = 1 n el Ω e where n el is the total number of elements in the mesh. We apply a multiscale overlapping decomposition to the velocity field and the weighting functions only in the narrow band Ω along the Dirichlet boundary. v = v + v (14) w = w + w (15) where v and v are the coarse and fine scale velocity fields, respectively. ŵ and w are the corresponding weighting functions for the coarse and fine scales, respectively. The coarse scale is associated with finite element spaces and the fine scale is represented by piecewise polynomials of sufficiently high order. Substituting Eqs (14) and (15) into the mixed weak form in Eqs (12) and (13) and using the linearity of the bilinear forms, we obtain two variational forms: The coarse-scale sub-problem governs the computable scales along Γ g , while the fine-scale sub-problem governs the residual-based error part along Γ g . Because we apply the scale decomposition only at Ω, the fine-scale problem Eq (18) is localized around the Dirichlet boundary Γ g . Accordingly, the fine scales v are assumed to be non-zero only within elements at the boundary Γ g and they vanish in the interior of the domain Ω.

Derivation of the fine-scale models
We assume that the fine-scale field is operational over the first layer of elements that are adjacent to the boundary Γ g . Consequently, it is nonzero on Ω, attaining maximum value at Γ g and becoming zero at distance ε = ℎ perpendicular to Γ g , where ℎ is the width of the element normal to the boundary. This modeling assumption allows the continuity of the overlapping coarse and fine fields. It also yields a simpler method for numerical integration and therefore for the approximation of the conditions at the edges. We approximate the finescale fields using edge functions b e ξ defined in the natural coordinate system ξ = ξ, η, ζ .
They are non-zero at the segment of Dirichlet boundaries Γ g e and zero at the other element boundaries. An example of edge function is shown in Figure 2 and Table 1 lists the edge functions for 2D and 3D elements. The fine scales are represented as follows: v = b e α (19) w = b e β (20) where α and β is the coefficient for the trial solutions and weighting functions of the fine scale velocities, respectively. We make a simplifying assumption that v is quasi-static so that the effect of its time-derivative is negligible, i.e., ∂v/ ∂t ≈ 0.
We linearize the fine-scale problem with respect to v by applying the variational derivative, . We gather the coarse-scale terms on the right-hand side and apply integration-by-parts and the divergence theorem to them. The resulting linearized problem for the incremental fine scales Δv is where r is the residual of the Euler Lagrange equations of the coarse scale formulation, and is given as The variation of the viscosity field is expressed as By substituting expressions Eqs (19), (20) and (23) into Eq (21), the fine-scale problem can be segregated into local problems over the elements Ω e that lie along the Dirichlet boundaries Γ g . b e β, ρ b e Δα ⋅ ∇v Ω e + b e β, ρv ⋅ ∇ b e Δα Ω e + ∇b e β, 2η γ˙ε b e Δα Ω e + ∇b e β, 4η′ γ˙γ˙− 1 ε v : ε b e Δα ε v Ω e = − b e β, r Ω e − b e β, σ ⋅ n + λ Γ g e − b e β, σ ⋅ n − h Γ ℎ e (24) We solve the linear problem for the fine-scale coefficients Δα and get the closed form expression for the fine-scale velocity field. Δv = − b e τ ∫ Ω e b e r dΩ + ∫ Γ g e b e σ ⋅ n + λ dΓ + ∫ Γ ℎ e b e σ ⋅ n − h dΓ (25) where τ = ∫ 4η′ γ˙γ˙− 1 ε v ⋅ ∇b e ⊗ ∇b e ⋅ ε v dΩ (26) The last boundary term in Eq (25) vanishes because the edge function b e is zero at the element edges on the traction boundary Γ ℎ e . To further simplify the fine-scale expression Eq (25), we employ three modeling assumptions [22,25]. (i) The edge bubble function b e is orthogonal to the coarse-scale residual r, so the domain-interior term in Eq (25) is neglected, i.e., ∫ Ω e b e r dΩ ≈ 0. Consequently, the fine-scale problem is driven by the boundary residual at the Dirichlet boundary Γ g e . (ii) The boundary residual is taken out of the integral by applying the mean-value theorem, i.e., ∫ Γ g e b e σ ⋅ n + λ dΓ ≈ σ ⋅ n + λ ∫ Γ g e b e dΓ . (iii) The edge function is averaged along the boundary, i.e., b e ≈ meas Γ g e −1 ∫ Γ g e b e dΓ where meas Γ g e is the length or the area of the element edge on Γ g e . With these modeling assumptions, we arrive at the fine-scale model at the boundary Γ g e that is driven by the boundary residual of the corresponding coarse scales.

Remark:
The last term in Eq (26) involves spatial gradient of nonlinear viscosity. This term would drop out for element-wise constant viscosity over a boundary element [19]. To keep the formulation general, we assume the viscosity to vary spatially over the element, and consequently we keep this term in the definition of the stability tensor.

Variational embedding in the coarse-scale formulation
In Section 3.2, we derived the fine-scale model for the incremental velocity field along the Dirichlet boundaries Γ g . The fine-scale velocity in Eq (27) is expressed in terms of the Lagrange multiplier λ which is still an unknown field. To derive a closed-form expression for λ, we substitute Eq (27) into Eq (17). Assuming a piecewise constant projection of λ along Γ g e , we solve Eq (29) locally to obtain the expression for the Lagrange multiplier λ. λ = − σ ⋅ n + τ g v − g (30) Notice that the Lagrange multiplier comprises the Cauchy traction and a penalty force to enforce the boundary condition. Substitution of Eq (30) in Eq (27) leads to an explicit form of the fine-scale model. (31) The fine scale in the current framework can be viewed as the residual or the error of the Dirichlet boundary condition Eq (3) at the boundary Γ g .
To embed the fine-scale model, we linearize the coarse-scale formulation Eq (16) with respect to the fine-scale velocity field v.
To convert the fine-scale terms that are integrated over the narrow band Ω into the boundary terms, we apply integration-by-parts. We keep only the boundary contributions and neglect the interior terms following the assumption that fine-scales are only active at the boundary Γ g . Substituting Eqs (30) and (31) in the linearized form Eq (32) yields the stabilized boundary formulation. The technical details of the steps for conversion from the domaininterior terms in Ω to the boundary terms on Γ g are provided in Appendix A. Furthermore, we add the residual-based interior stabilization χ, τr Ω for flows of the shear-rate dependent fluids developed in [18,19]. Hereon, we drop the superposed hats from the coarse-scale velocity field and its weighting function. Note that the stability tensor τ is computed using the regular bubble function b ξ , while the boundary stability tensor τ g is computed using the edge function b e ξ .
In Eq (33), the stabilized formulation is augmented by the weakly imposed boundary condition. The boundary terms integrated over Γ g weakly enforce the Dirichlet boundary condition. The first three boundary terms ensure the consistency of the formulation, which is required for optimal convergence. The fourth boundary term emanates from the linearization of the shear stress term of the fine-scale problem, which disappears for the Newtonian fluid case. The last boundary term stabilizes the boundary formulation, where the stability tensor is self-calculated without user-parameter.

Remark:
The derived formulation Eq (33) can be used with any shear-rate dependent fluid model. In this work, we have employed the power-law and the Carreau-Yasuda models. The formulation also accommodates the degenerate case of Newtonian fluids where the viscosity becomes constant and the viscosity-gradient term drops out.
Remark: An important point to note is that all the boundary terms are variationally derived and therefore the resulting formulation Eq (33) is variationally consistent.

Remark:
The fourth boundary term in the formulation Eq (33) involving the gradient of viscosity emanates from the linearization of the fine-scale problem. The effects of this term on the convergence of the formulation are discussed in Section 4.1.2. In the numerical experiment of shear-thinning flow in a curved tube discussed in Section 4.3.2, simulations without this term diverge unless the time-step size is substantially reduced. It shows that this term helps in preserving unconditional stability of the implicit time integration methods.

Remark:
The fourth boundary term in the formulation Eq (33) involves the inverse of the shear-rate, which can cause singularity at γ ≈ 0. The boundedness of this term is shown in Appendix B. In our numerical implementation, this boundary term is deactivated if γ < ε where ε = 10 −16 .

Numerical results
We have implemented the stabilized boundary formulation Eq (33) and its consistent tangent tensors using linear quadrilateral elements for 2D problems and quadratic tetrahedral elements for 3D problems. The generalized-α method is used for time integration. It is an implicit, second-order accurate scheme with a free parameter 0 ≤ ρ ∞ ≤ 1 that controls the damping of high frequency modes. In this work, we have used ρ ∞ = 0.5 for all the transient test cases. The nonlinear problems are solved using the Newton-Raphson method with the consistent tangent. We solve the linear system of equations using the direct solver for convergence tests and 2D cavity flow, and the GMRES solver with additive Schwarz preconditioner for the 3D arterial flow test cases. Note that the velocity field is divergence-free. The corresponding body force for the exact velocity and pressure fields is given in Appendix A in [18].
The density of the fluid is = 1.0 kg ⋅ m −3 . The material parameters for Carreau-Yasuda models are set as μ ∞ = 0.00345 Pa ⋅ s, μ 0 = 0.056 Pa ⋅ s, λ = 1.902, n = 0.22 and a = 1.25. Based on the maximum velocity and the viscosity at the infinite shear-rate μ ∞ , the Reynolds number is Re = 435. The flow is driven by the body force given in [18]. The values of the exact velocities are applied as nodal coefficients on all boundaries for the case of strongly imposed boundary conditions, while they are enforced via the formulation in Eq (33) for the weakly imposed boundary conditions. Figure 3 presents the velocity and pressure fields computed on the finest mesh with the weakly imposed boundary conditions. Figure 4 shows the convergence rates for the L 2 -norm and the H 1 seminorm of error that shows optimal convergence rates for the velocity field and its divergence. Both strong and weak boundary conditions show equivalent convergence properties for the finer meshes, while the weak boundary condition results in smaller error for the coarse mesh. This shows that the weakly imposed boundary conditions help improve the accuracy of calculations on cruder spatial discretizations.

Domain with penetration boundaries-We now develop a harder test case for
the convergence-rate study presented in Section 4.1.1. We alter the computational domain by cutting the original domain right through the vortex such that −0.5 ≤ x, y, z ≤ 0.25. Since exact solution for this problem exists, which is given in Eq (36), the boundary conditions at this bounding surface can be applied both strongly as well as weakly. It is important to realize that in this problem where the domain boundary cuts through the vortex structure as shown in Figure 5, imposing the boundary conditions becomes more involved. Figure 6(a) compares the convergence rates for the strongly and weakly imposed boundary conditions. We observe that for the case of strongly imposed boundary conditions the L 2 -norm of error for the pressure field does not converge for the finest mesh, while the errors for the weakly imposed boundary conditions uniformly converge for all the meshes. Figure 6(b) compares the convergence properties between the complete formulation Eq (33) and the incomplete formulation that is obtained by deactivating the boundary term that accounts for the viscosity gradient, i.e., the fourth boundary term in Eq (33). We realize that not including the boundary term with the shear-rate effects leads to catastrophic divergence of the method for the finest mesh, while keeping this term in the formulation leads to uniform convergence all through the mesh refinement. This numerical study justifies the crucial importance of this term in practical applications with non-Newtonian shear-rate fluids.

Lid-driven cavity flow
Steady lid-driven cavity flow is chosen to verify the spatial accuracy of the formulation presented in Eq (33) for both Newtonian and non-Newtonian fluids. We compare the results with the benchmark data [15,[26][27][28] and with the numerical results for the strongly imposed boundary conditions. Figure 7 shows the biunit domain with the Dirichlet boundary conditions for the velocity field applied along all the faces. A constant unit velocity U = 1 is applied in the x-direction at the top surface while no-slip boundary condition is applied at the other three surfaces. The zero pressure condition is specified at the left bottom corner to eliminate the arbitrary constant. We test both the Newtonian and the shear-thinning models as presented below.

Newtonian fluids-For
Newtonian fluids, Reynolds number is defined based on the lid-velocity U and domain width L. We consider three cases of Reynolds number Re = 1000,5000 by setting the density ρ = 1 and the viscosity μ = 10 −3 , 2 × 10 −4 , respectively. The computational domain is discretized using uniformly distributed N × N (where N = 40,80) quadrilateral elements. Figure 8 presents the profiles of the horizontal velocity along the vertical center line at x = 0.5. We compare the computed results of strongly and weakly imposed boundary conditions with the reference results that are computed on the stretched meshes [26,26]. We observe that the weakly imposed boundary conditions successfully capture the reference profile and they outperform the strongly imposed boundary conditions on all the meshes.

Remark:
The power-law model gives rise to singularity when γ ≈ 0. To keep the viscosity bounded at very low shear-rates, we impose a lower cap on the shear-rate, γ ≥ ε, where we have employed ε = 10 −16 . This treatment is equivalent to applying an upper bound for the viscosity, η γ ≤ η ε . For the shear-thinning power-law fluids in Section 4.2.2, the viscosity is bounded by η γ ≤ 200,000.
Figures 9-12 present the comparison of the computed velocity profiles on the stretched and on the uniform meshes with the reference results [15,28]. The reference results for the power-law and the Carreau-Yasuda models are computed on 180 × 180 and 128 × 128 meshes, respectively. We achieve good comparative results on meshes of different resolution for the power-law and Carreau-Yasuda models which indicates that the proposed method can handle different non-Newtonian constitutive equations. For the case of stretched meshes, all mesh configurations yield very accurate results in the boundary layer region, as shown in Figures 9 and 11. For the case of uniform meshes, the numerical results for both strongly and weakly imposed boundary conditions approach the reference results as a function of mesh refinement, as shown in Figures 10 and 12. However, the weakly imposed boundary conditions result in higher accuracy in the boundary layer region as compared to the strongly imposed boundary conditions. This feature demonstrates the enhanced mathematical attributes of the weakly imposed boundary conditions on the uniform meshes that may not have optimal mesh resolution for the boundary-layer flows. Figure 13 shows the magnitude of the non-dimensionalized stabilization tensor ∥ τ g ∥ℎ/μ ∞ computed on the uniform meshes along the left vertical boundary for the case of Carreau-Yasuda model. The stabilization tensor τ g defined in Eq (28) is calculated over the layer of elements adjacent to the boundary Γ g . The magnitude of the tensor is computed as ∥ τ g ∥ = τ g : τ g , ℎ is the length of the element edge, and μ ∞ is the asymptotic viscosity at the infinite shear-rate in Carreau-Yasuda model. As the meshes are refined, the width of the spatial domain comprising the first layer of elements adjacent to Γ g reduces, thereby asymptoting to the dimension n sd − 1. Figure 13 shows that all the meshes yield similar spatial distribution of the non-dimensional stability parameter and they gradually converge to the results from the finest mesh that yields a better representation of the notion of the edge.

Flow in a curved tube
This test case is a model for flow in an idealized femoral artery [29] where we investigate both Newtonian and non-Newtonian constitutive equations. The geometry of the curved tube and its dimensions are shown in Figure 14. The diameter of the tube is D = 0.8 cm, the length of the straight segment is L = 1.6 cm, the radius of curvature is R = 2.4 cm, and the angle from the inlet is θ = 0 90 ∘ . The parabolic inflow velocity profile is prescribed at the bottom surface. The no-slip boundary condition is strongly or weakly imposed over the cylindrical surface, and the traction-free condition is applied at the outflow surface. We employ two successively refined quadratic tetrahedral meshes where the refined mesh is generated by reducing the element size by half. The description of the two meshes is given in Table 2.

4.3.1.
Verification of the method with the Newtonian fluid model-We first test our method with Newtonian fluids and compare the computed results with the experimental data [29] that was obtained via laser-Doppler velocity measurements. The density of fluid is ρ = 0.8028 g/cm 3 and the dynamic viscosity is μ = 0.1 dyne ⋅s/cm 2 . The average velocity at the inlet is U = 109 cm/s, which yields Reynolds number Re = ρUD/μ = 700. Figure 15 presents the streamwise velocity along the cross sections subtended at various angles as shown in Figure 14. A good agreement between the computed results and experimental data is attained on the fine mesh. The coarse mesh has only eight quadratic elements along the diameter, with no-slip condition at the ends of the diameter, thereby giving rise to boundary layers at either end. Despite being relatively crude, the mesh is able to capture the profile quite accurately and this is attributed to the fine-scale modeling feature of the method. We also observe that both strongly and weakly imposed boundary conditions produce similar results on each mesh.

Shear-thinning non-Newtonian blood flow-
In this section, we consider pulsatility of blood flow in addition to the shear-thinning effects of blood. We employ the Carreau-Yasuda model in Eq (7) with the coefficients given in Table 3 that are taken from [13]. The density of blood is ρ = 1.06 g/cm 3 . The prescribed timevarying velocity at the center of the inlet is presented in Figure 16. We simulate three cardiac cycles with a time step size of Δt = 0.01 s and take the numerical results from the last cycle for comparative study. Figure 17 presents the computed velocity field on the fine mesh by weakly imposing the no-slip condition at the wall. It shows that the boundary terms in Eq (33) accurately enforce the no-slip condition that leads to the boundary layer near the wall. Figure 18 compares the computed streamwise velocity at various angles between the strongly and weakly imposed boundary conditions. Figure 19 shows the wall shear stress (WSS) along the outermost curve of the bent tube. The WSS is averaged in time for a cardiac cycle. The comparisons of the velocity and WSS verify that both methods for imposing the essential boundary condition produce comparable results on the coarse and fine meshes. Table 4 summarizes the range of eigenvalues of the stabilization tensor τ for the case of weakly imposed boundary conditions. Since the eigenvalues remain positive, it shows that the stability tensor also stays bounded and positive-definite.
A unique attribute of the proposed method is the variationally derived stabilization tensor τ g that enforces the Dirichlet boundary condition. It is explicitly defined in Eq (28) and it results in element-wise optimal values for τ g that are automatically calculated locally. Figure  20 displays the spatial distribution of the magnitude of the stabilization tensor τ g computed at two time points T A and T B in a typical cardiac cycle. It shows that the value of τ g tends to be larger where the velocity and its gradient are higher on the top surface and near the outlet. It also shows that the value of τ g varies with time and is larger at the time point T A than at the time point T B because the inflow velocity is also higher at T A .

Blood flow in a representative patient-specific arterial model
In this test case we take the method for weakly imposed boundary conditions to patientspecific aortic and femoral arteries of clinical relevance. We constructed a patient-specific geometric model using the cross-sectional CT scan images and the 3D geometry of the arterial tree branching from the thoracic aorta to the femoral arteries as shown in Figure 21.
The geometric model was discretized using quadratic tetrahedral elements where the number of nodes and elements in the mesh are 586,116 and 388,810, respectively. We simulated three cardiac cycles with a time increment of Δt = 0.005 s, which corresponds to 220 time steps during a typical cardiac cycle.
The non-Newtonian effects of blood are accounted for via the Carreau-Yasuda model with parameters given in Table 3. The density of blood is ρ = 1.06 g/cm 3 . The flowrate shown in Figure 22 is specified at the inflow surface of the thoracic artery by imposing the parabolic distribution of velocity. The artery wall is considered as no-slip surface. The resistance boundary condition in Eq (38) that produces physiological pressure waveform at the outlets is imposed at all the outlets to incorporate the downstream resistive effects of the arteries [17,30].
where R is the resistance parameter and p 0 = 80 mmHg is the constant downstream pressure.
The resistance parameter for each branch is tuned based on the flow-rate in the branches, which is calculated with the traction-free conditions, to obtain a pressure amplitude between 80-110 mmHg at the outlets. Table 5 summarizes the value of the resistance parameter for each branch.
Non-physical backflow at the outlet is often encountered in cardiovascular simulations, and it can destabilize the solution. To address this issue, we apply backflow stabilization [31,32] to the outflow boundary surface by adding the following term to the left-hand side of Eq (33).
This term is only activated where the velocity is pointed inwards to the fluid domain at the outflow surface.

Remark:
In this patient-specific test case, where no-slip condition is weakly imposed all over the arterial wall surface, the derived stabilization tensor τ g defined in Eq (28) is multiplied by a constant α = 20, which leads to better convergence in the non-linear solution. Figure 23 compares the flowrate at the outlets of some representative branches for the cases of strongly and weakly imposed essential boundary conditions, and a good agreement between the computed solutions is observed. Figure 24 shows that the computed pressure profiles with the strong and weak imposition of boundary conditions closely match for the inflow boundary and the outflow boundary, respectively. This test case verifies that the weakly imposed boundary condition produces equivalent blood flow patterns in the arterial system as are produced by the strongly imposed boundary conditions. Figures 25-27 present the numerical results for blood flow in the patient-specific arteries computed via the weakly imposed boundary conditions. Figure 25 shows the volume rendering of the magnitude of velocity field. Figure 26 shows the instantaneous snapshots of the viscosity field at the peak systole and the mid diastole, where the viscosity tends to be lower at the systole due to higher blood velocity. We observe the wide spread of the high-viscosity region during the diastole, especially around the thoracic aorta near the inlet, and the bifurcating parts of arteries. The high viscosity shows the propensity of blood to coagulate. Figure 27 shows the spatial distribution of the wall shear stress (WSS), which is the tangential component of the stress acting on the arterial wall, and is considered a significant factor for the progression of arterial disease. We want to note that the viscosity and WSS data is difficult to obtain via in vivo experiments. Computational methods can enable us to identify the regions where viscosity or WSS are higher, and therefore can serve as virtual platforms for patient specific care and surgical planning.

Conclusions
We have presented a stabilized method for weakly enforcing the Dirichlet boundary conditions for non-Newtonian shear-rate dependent fluids. The consistency and stabilization terms are derived by locally resolving the fine-scale variational equation along the Dirichlet boundary. For shear-rate fluids, additional boundary terms appear that are not present in the Nitsche type approaches that are otherwise employed for weak enforcement of Dirichlet boundary conditions. One of these terms is the viscosity gradient term and is a function of the shear-rate, while the other term is a zeroth-order term. The structure of the edge stabilization tensor that weights the boundary terms appears naturally, and it is free of user defined parameters. Employing edge functions the edge stabilization tensor is calculated automatically, and it adaptively adjusts itself to the magnitude of the boundary residual. Another significant advantage of the method for weakly imposed boundary conditions originates from the flexibility of not requiring nodally matched discretizations at bloodartery interaction surfaces to enforce continuity of primary fields in blood-tissue interaction problems. This relaxation in the generation of meshes for the blood and tissue sub-regions will be beneficial in complex geometries arising in patient-specific models.
The convergence rate tests show that the weakly imposed boundary conditions achieve optimal order of convergence for the velocity and pressure fields in the norms considered. Furthermore, the boundary term containing viscosity-derivative and shear-rate of the fluid is essential to obtain stable and convergent solution as a function of mesh refinement. The weakly imposed boundary conditions for 2D cavity flow outperform the strongly imposed boundary conditions in resolving the boundary layers and yield higher spatial accuracy of velocity gradients. The test case of 3D bent tube verifies the spatial accuracy of the solution in comparison with the experimental data for Newtonian fluids, and via comparison with the case of strongly imposed boundary conditions for non-Newtonian fluids. The distribution of the magnitude of the stability parameter at the boundary highlights that it adaptively adjusts itself spatially and temporally in response to the magnitude of the residual at the boundary. The test case of patient-specific arteries combines the various mathematical attributes of the method and shows its potential for the clinically relevant applications.
The following example shows the application of these transformations for the viscous term. Assuming that the incremental fine-scale field is active, Δv ≠ 0, only at the Dirichlet boundaries Γ g (Assumption 1), the boundary term on Γ is relegated to the boundary term on Γ g . Additionally, we neglect the domain-interior terms including Δv based on the assumption that the fine-scale field is vanishingly small, Δv ≈ 0, in the domain interior near the boundary Ω (Assumption 2). B = ∫ Γ g 2η γ˙ε w ⋅ n ⋅ Δv dΓ  Profiles of the streamwise velocity across the cross-sections at various angles (The radial distance is normalized with respect to the radius of the tube, and the velocity is normalized with respect to the average velocity at the inlet).  The velocity waveform applied at the center of the inflow surface (Time points T A and T B correspond to the maximum and minimum inflow velocity). Numerical results for the velocity field on the fine mesh by weakly imposed no-slip condition. Profiles of the streamwise velocity across the cross-sections at various angles at the time point T A where the inflow velocity is the maximum (The radial distance is normalized with respect to the radius of the tube).  Time-averaged wall shear stress (WSS) for a cardiac cycle along the outermost circular curve in the bent tube. Kang   Spatial distribution of the magnitude of the non-dimensionalized stabilization tensor ∥ τ g ∥ℎ/μ ∞ computed on the intermediate mesh at two time points T A and T B indicated in Figure 16. The magnitude of the tensor is computed as ∥ τ g ∥ = τ g : τ g , ℎ is the element length-scale, and μ ∞ is the asymptotic viscosity at the infinite shear-rate in Carreau-Yasuda model.

Author Manuscript
Author Manuscript Comparison of the pressure at the inflow and outflow surfaces between the strongly and weakly imposed BC cases. (The outflow pressure is computed at outlet 7 marked in Figure  21). Kang    Viscosity contours computed via the weakly imposed boundary condition at the peak systole and the mid diastole. (Unit: dyne · s/cm 2 ). Wall shear stress (WSS) computed via the weakly imposed boundary condition at the peak systole and the mid diastole. (Unit: dyne · s/cm 2 ). Kang  Description of the meshes for the curved tube (ℎ is the side length of the element).