Artificial intelligence in acupuncture: A bibliometric study

1.   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 blood-artery 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 ^[16]. 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.

2.   Governing equations and non-Newtonian fluid models

2.1.   The strong form of governing equations

Let $\Omega \subset {\mathbb{R}^{{n_{{\kern 1pt} {\rm{sd}}}}}}$ be an open bounded region with piecewise smooth boundary $\Gamma$. The number of spatial dimensions, ${n_{{\kern 1pt} {\rm{sd}}}}$, is equal to 2 or 3. The domain boundary assumes the usual split $\Gamma = {\Gamma _g} \cup {\Gamma _h}$ and ${\Gamma _g} \cap {\Gamma _h} = \emptyset$, where ${\Gamma _g}$ and ${\Gamma _h}$ are parts of the boundary with essential and natural boundary conditions, respectively. The governing equations and boundary conditions are given as:

where ${\boldsymbol{v}}$ and $p$ are the fluid velocity and pressure, respectively. $\rho$ is the density of the fluid, ${{\boldsymbol{\sigma }}_v}$ is the deviatoric stress tensor, $\mathit{\boldsymbol{\sigma }}$ is the Cauchy stress tensor defined as $\mathit{\boldsymbol{\sigma }} = - p{\boldsymbol{I}} + {\mathit{\boldsymbol{\sigma }}_v}$, ${\boldsymbol{I}}$ is the identity tensor, ${\boldsymbol{f}}$ is the body force per unit mass, ${\boldsymbol{g}}$ is the prescribed velocity on the boundary ${\Gamma _g}$, ${\boldsymbol{h}}$ is the prescribed traction on the boundary ${\Gamma _h}$, and ${\boldsymbol{n}}$ is the unit outward normal to the boundary of the domain $\Omega$. 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 ${\boldsymbol{\varepsilon }} ({\boldsymbol{v}})$ is the rate-of-deformation tensor defined as $\mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}}) = \frac{1}{2}({\kern 1pt} {\kern 1pt} \nabla {\boldsymbol{v}} + \nabla {{\boldsymbol{v}}^T})$, $\dot \gamma$ is the shear-rate defined as $\dot \gamma = \sqrt {2\mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}}):\mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}})}$, and $\eta ({\kern 1pt} {\kern 1pt} \dot \gamma {\kern 1pt} {\kern 1pt})$ is the effective viscosity which is a function of the shear-rate.

2.2.   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:

where $\mu$ 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

where ${\mu _{{\kern 1pt} 0}}$ and ${\mu _{{\kern 1pt} \infty }}$ are asymptotic viscosities at zero and infinite shear-rate, respectively, and $a$, $n$ and $\lambda$ 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 ${\mu _{{\kern 1pt} 0}} = {\mu _{{\kern 1pt} \infty }}$.

2.3.   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 ${\Gamma _g}$. Let ${\boldsymbol{w}}$ and $q$ denote the weighting functions for the velocity and the pressure fields, and $\mathit{\boldsymbol{\psi }}$ 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 $\left\{ {{\kern 1pt} {\boldsymbol{v}}, p{\kern 1pt} } \right\} \in {\mathcal{S}} \times {\mathcal{P}}$, $\mathit{\boldsymbol{\lambda }} \in {\mathcal{Q}}$ such that for all $\left\{ {{\kern 1pt} {\boldsymbol{w}}, q{\kern 1pt} } \right\} \in {\mathcal{V}} \times {\mathcal{P}}$, $\mathit{\boldsymbol{\psi }} \in {\mathcal{Q}}$:

where ${({\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt}, {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdot {\kern 1pt} {\kern 1pt})_\Omega } = \int_\Omega {(\, {\kern 1pt} \cdot \, {\kern 1pt}){\kern 1pt} {\kern 1pt} d\Omega }$ is the ${L_{{\kern 1pt} {\rm{2}}}}{\rm{(}}{\kern 1pt} {\kern 1pt} \Omega {\kern 1pt} {\kern 1pt} {\rm{)}}$ 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 $\mathit{\boldsymbol{\lambda }}$ which acts as the numerical flux or the traction at ${\Gamma _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 $\mathit{\boldsymbol{\lambda }}$ that will eliminate the auxiliary unknown field from the formulation.

3.   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 $\tilde \Omega \subset \Omega$ along the Dirichlet boundary ${\Gamma _g}$ shown in Figure 1. Overall procedure underlying the derivation comprises three steps, (ⅰ) split the problem into coarse- and fine-scale sub-problems, (ⅱ) solve the fine-scale problem along the boundary ${\Gamma _g}$, and (ⅲ) embed the fine-scale models into the coarse-scale formulation.

3.1.   Multiscale decomposition

We discretize the domain $\Omega$ into disjoint elements ${\Omega ^e}$ with element boundaries ${\Gamma ^{{\kern 1pt} e}}$, such that $\Omega \, \, = \, \, \cup {\kern 1pt} _{e{\kern 1pt} {\kern 1pt} = {\kern 1pt} {\kern 1pt} 1}^{{n_{{\kern 1pt} {\rm{e}}{\kern 1pt} {\rm{l}}}}}\, {\kern 1pt} {\Omega ^e}$ where ${n_{{\kern 1pt} {\rm{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 $\tilde \Omega$ along the Dirichlet boundary.

where ${\boldsymbol{\hat v}}$ and ${\boldsymbol{\tilde v}}$ are the coarse and fine scale velocity fields, respectively. ${\boldsymbol{\hat w}}$ and ${\boldsymbol{\tilde 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:

Coarse-scale sub-problem

Fine-scale sub-problem

The coarse-scale sub-problem governs the computable scales along ${\Gamma _g}$, while the fine-scale sub-problem governs the residual-based error part along ${\Gamma _g}$. Because we apply the scale decomposition only at $\tilde \Omega$, the fine-scale problem Eq (18) is localized around the Dirichlet boundary ${\Gamma _g}$. Accordingly, the fine scales ${\boldsymbol{\tilde v}}$ are assumed to be non-zero only within elements at the boundary ${\Gamma _g}$ and they vanish in the interior of the domain $\Omega$.

3.2.   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 ${\Gamma _g}$. Consequently, it is nonzero on $\tilde \Omega$, attaining maximum value at ${\Gamma _g}$ and becoming zero at distance $\varepsilon = h$ perpendicular to ${\Gamma _g}$, where $h$ 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 fine-scale fields using edge functions ${b^e}({\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{\xi }}{\kern 1pt} {\kern 1pt} {\kern 1pt})$ defined in the natural coordinate system $\mathit{\boldsymbol{\xi }} = ({\kern 1pt} {\kern 1pt} \xi, \; \eta, \; \zeta {\kern 1pt} {\kern 1pt})$. They are non-zero at the segment of Dirichlet boundaries $\Gamma _g^{{\kern 1pt} 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:

where $\mathit{\boldsymbol{\alpha }}$ and $\mathit{\boldsymbol{\beta }}$ is the coefficient for the trial solutions and weighting functions of the fine scale velocities, respectively. We make a simplifying assumption that ${\boldsymbol{\tilde v}}$ is quasi-static so that the effect of its time-derivative is negligible, i.e., $\partial {\kern 1pt} {\boldsymbol{\tilde v}}/\partial {\kern 1pt} t \approx {\bf{0}}$.

We linearize the fine-scale problem with respect to ${\boldsymbol{\tilde v}}$ by applying the variational derivative, ${D_{\boldsymbol{v}}}[{\kern 1pt} {\kern 1pt} {\boldsymbol{G}}{\kern 1pt} {\kern 1pt} {\kern 1pt}] \cdot \Delta {\boldsymbol{v}}\; \, = \; \, \frac{\partial }{{\partial \varepsilon }}{\left[{{\kern 1pt} {\boldsymbol{G}}{\kern 1pt} {\kern 1pt} (\, {\boldsymbol{v}}\; + \; \varepsilon {\kern 1pt} \Delta {\boldsymbol{v}}\, ){\kern 1pt} } \right]_{{\kern 1pt} \varepsilon {\kern 1pt} = {\kern 1pt} {\kern 1pt} 0}}$. 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 $\Delta {\boldsymbol{\tilde v}}$ is

where ${\boldsymbol{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 ${\Omega ^e}$ that lie along the Dirichlet boundaries ${\Gamma _g}$.

We solve the linear problem for the fine-scale coefficients $\Delta \mathit{\boldsymbol{\alpha }}$ and get the closed form expression for the fine-scale velocity field.

where

The last boundary term in Eq (25) vanishes because the edge function ${b^e}$ is zero at the element edges on the traction boundary $\Gamma _h^{{\kern 1pt} e}$. To further simplify the fine-scale expression Eq (25), we employ three modeling assumptions ^[22,25]. (ⅰ) The edge bubble function ${b^e}$ is orthogonal to the coarse-scale residual ${\boldsymbol{r}}$, so the domain-interior term in Eq (25) is neglected, i.e., $\int_{{\Omega ^e}} {{b^e}{\boldsymbol{r}}\, d\Omega } \approx 0$. Consequently, the fine-scale problem is driven by the boundary residual at the Dirichlet boundary $\Gamma _g^{{\kern 1pt} e}$. (ⅱ) The boundary residual is taken out of the integral by applying the mean-value theorem, i.e., $\int_{\Gamma _g^{{\kern 1pt} e}} {{b^e}(\mathit{\boldsymbol{\sigma }} \cdot {\boldsymbol{n}} + {\boldsymbol{\lambda }} \, ){\kern 1pt} {\kern 1pt} d\Gamma } \approx (\mathit{\boldsymbol{\sigma }} \cdot {\boldsymbol{n}} + {\boldsymbol{\lambda }} \, )\int_{\Gamma _g^{{\kern 1pt} e}} {{b^e}{\kern 1pt} d\Gamma }$. (ⅲ) The edge function is averaged along the boundary, i.e., ${b^e} \approx {{\rm{[meas}}\, (\, \Gamma _g^{{\kern 1pt} e})]^{{\kern 1pt} - {\kern 1pt} {\kern 1pt} 1}}\int_{\Gamma _g^{{\kern 1pt} e}} {{b^e}d\Gamma }$ where ${\rm{meas}}\, (\, \Gamma _g^{{\kern 1pt} e})$ is the length or the area of the element edge on $\Gamma _g^{{\kern 1pt} e}$. With these modeling assumptions, we arrive at the fine-scale model at the boundary $\Gamma _g^{{\kern 1pt} e}$ that is driven by the boundary residual of the corresponding coarse scales.

where the stability tensor is defined as

where $\tilde \tau$ is given by Eq (26).

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 ^[20]. 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.

3.3.   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 ${\Gamma _g}$. The fine-scale velocity in Eq (27) is expressed in terms of the Lagrange multiplier $\mathit{\boldsymbol{\lambda }}$ which is still an unknown field. To derive a closed-form expression for $\mathit{\boldsymbol{\lambda }}$, we substitute Eq (27) into Eq (17).

Assuming a piecewise constant projection of $\mathit{\boldsymbol{\lambda }}$ along $\Gamma _g^{{\kern 1pt} e}$, we solve Eq (29) locally to obtain the expression for the Lagrange multiplier $\mathit{\boldsymbol{\lambda }}$.

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.

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 ${\Gamma _g}$.

To embed the fine-scale model, we linearize the coarse-scale formulation Eq (16) with respect to the fine-scale velocity field ${\boldsymbol{\tilde v}}$.

To convert the fine-scale terms that are integrated over the narrow band $\tilde \Omega$ 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 ${\Gamma _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 domain-interior terms in $\tilde \Omega$ to the boundary terms on ${\Gamma _g}$ are provided in Appendix A. Furthermore, we add the residual-based interior stabilization ${\left({\mathit{\boldsymbol{\chi }}{\kern 1pt}, \mathit{\boldsymbol{\tau }}{\kern 1pt} {\boldsymbol{r}}} \right)_\Omega }$ for flows of the shear-rate dependent fluids developed in ^[18,20]. Hereon, we drop the superposed hats from the coarse-scale velocity field and its weighting function.

where

Note that the stability tensor $\mathit{\boldsymbol{\tau }}$ is computed using the regular bubble function $b(\mathit{\boldsymbol{\xi }}\, )$, while the boundary stability tensor ${\mathit{\boldsymbol{\tau }}_g}$ is computed using the edge function ${b^e}(\mathit{\boldsymbol{\xi }}\, )$.

In Eq (33), the stabilized formulation is augmented by the weakly imposed boundary condition. The boundary terms integrated over ${\Gamma _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 $\dot \gamma \approx 0$. The boundedness of this term is shown in Appendix B. In our numerical implementation, this boundary term is deactivated if $\dot \gamma < \varepsilon$ where $\varepsilon = {10^{ - 16}}$.

4.   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\; {\kern 1pt} \leqslant \; {\kern 1pt} {\rho _\infty } \leqslant \; {\kern 1pt} 1$ that controls the damping of high frequency modes. In this work, we have used ${\rho _\infty } = 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.

4.1.   Rate of convergence study

4.1.1.   Domain with no-penetration boundaries

This test case ^[18] investigates the convergence-rates for the strongly and weakly imposed boundary conditions. The domain is a 3D block, $- {\kern 1pt} {\kern 1pt} 0.5\; {\kern 1pt} \leqslant \; {\kern 1pt} x, y, z\; {\kern 1pt} \leqslant \; {\kern 1pt} 0.5$, and it is discretized using evenly spaced quadratic tetrahedral elements. The number of elements per edge is $N =$2, 4, 8, 16, 32. The exact solution for the velocity and pressure fields is given as follows

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 $\rho = 1.0{\rm{ kg}} \cdot {{\rm{m}}^{ - {\kern 1pt} {\kern 1pt} 3}}$. The material parameters for Carreau-Yasuda models are set as ${\mu _{{\kern 1pt} \infty }} = 0.00345{\rm{ Pa}} \cdot {\rm{s}}$, ${\mu _{{\kern 1pt} 0}} = 0.056{\rm{ Pa}} \cdot {\rm{s}}$, $\lambda = 1.902$, $n = 0.22$ and $a = 1.25$. Based on the maximum velocity and the viscosity at the infinite shear-rate ${\mu _{{\kern 1pt} \infty }}$, 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_{{\kern 1pt} 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.

4.1.2.   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 $- {\kern 1pt} {\kern 1pt} 0.5\; {\kern 1pt} \leqslant \; {\kern 1pt} x, y, z\; {\kern 1pt} \leqslant \; {\kern 1pt} 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_{{\kern 1pt} 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.

4.2.   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 ^{[16,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.

4.2.1.   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, {\rm{ }}5000$ by setting the density $\rho = 1$ and the viscosity $\mu = {10^{ - {\kern 1pt} 3}}$, ${\rm{2}} \times {10^{ - {\kern 1pt} 4}}$, respectively. The computational domain is discretized using uniformly distributed $N \times N$ (where $N = {\rm{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,27]. 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.

4.2.2.   Shear-thinning non-Newtonian fluids

In this section, we employ the power-law and the Carreau-Yasuda models given in Eqs (6) and (7), respectively, for the shear-thinning fluids. The Reynolds number for the power-law fluids is defined by $Re = \rho \, {U^{2 - n}}{L^{{\kern 1pt} n}}/\mu$, as given in reference ^[16]. We set the exponent $n = 0.5$, the density $\rho = 1$, and the viscosity parameter $\mu = 2 \times {10^{ - {\kern 1pt} 3}}$, which yields $Re = 500$. The material parameters for Carreau-Yasuda fluids are set equal to ${\mu _{{\kern 1pt} \infty }} = 0.00345$, ${\mu _{{\kern 1pt} 0}} = 0.056$, $\lambda = 1.902$, $n = 0.22$, $a = 1.25$, and $\rho = 1$, which yields the Reynolds number $Re = 290$ based on ${\mu _{{\kern 1pt} \infty }}$. We employ (i) stretched, and (ii) uniform quadrilateral meshes, where the number of elements along an edge is $N =$40, 80,160 for the power-law model, and $N =$20, 40, 80 for the Carreau-Yasuda model.

Remark: The power-law model gives rise to singularity when $\dot \gamma \approx 0$. To keep the viscosity bounded at very low shear-rates, we impose a lower cap on the shear-rate, $\dot \gamma \geqslant \varepsilon$, where we have employed $\varepsilon = {10^{ - 16}}$. This treatment is equivalent to applying an upper bound for the viscosity, $\eta (\, \dot \gamma \, ) \leqslant \eta (\varepsilon)$. For the shear-thinning power-law fluids in Section 4.2.2, the viscosity is bounded by $\eta (\, \dot \gamma \, ) \leqslant 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 ^[16,28]. The reference results for the power-law and the Carreau-Yasuda models are computed on $180 \times 180$ and $128 \times 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 $\left\| {{\mathit{\boldsymbol{\tau }}_g}} \right\|h/{\mu _{{\kern 1pt} \infty }}$ computed on the uniform meshes along the left vertical boundary for the case of Carreau-Yasuda model. The stabilization tensor ${\tau _g}$ defined in Eq (28) is calculated over the layer of elements adjacent to the boundary ${\Gamma _g}$. The magnitude of the tensor is computed as $\left\| {{\mathit{\boldsymbol{\tau }}_g}} \right\| = \sqrt {{\mathit{\boldsymbol{\tau }}_g}:{\mathit{\boldsymbol{\tau }}_g}}$, $h$ is the length of the element edge, and ${\mu _{{\kern 1pt} \infty }}$ 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 ${\Gamma _g}$ reduces, thereby asymptoting to the dimension ${n_{{\kern 1pt} {\rm{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.

4.3.   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{\rm{ cm}}$, the length of the straight segment is $L = 1.6{\rm{ cm}}$, the radius of curvature is $R = 2.4{\rm{ cm}}$, and the angle from the inlet is $\theta = 0 \sim 90^\circ$. 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 $\rho = 0.8028{\rm{ g}}{\kern 1pt} {\rm{/c}}{{\rm{m}}^3}$ and the dynamic viscosity is $\mu = 0.1{\rm{ dyne}} \cdot {\rm{s}}{\kern 1pt} {\rm{/c}}{{\rm{m}}^2}$. The average velocity at the inlet is $U = 109{\rm{ cm}}{\kern 1pt} {\rm{/s}}$, which yields Reynolds number $Re = \rho {\kern 1pt} UD/\mu = 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.

4.3.2.   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 $\rho = 1.06{\rm{ g}}{\kern 1pt} {\rm{/c}}{{\rm{m}}^3}$. The prescribed time-varying velocity at the center of the inlet is presented in Figure 16. We simulate three cardiac cycles with a time step size of $\Delta t = 0.01{\rm{ 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 $\tau$ 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 ${\mathit{\boldsymbol{\tau }}_g}$ that enforces the Dirichlet boundary condition. It is explicitly defined in Eq (28) and it results in element-wise optimal values for ${\mathit{\boldsymbol{\tau }}_g}$ that are automatically calculated locally. Figure 20 displays the spatial distribution of the magnitude of the stabilization tensor ${\mathit{\boldsymbol{\tau }}_g}$ computed at two time points ${T_A}$ and ${T_B}$ in a typical cardiac cycle. It shows that the value of ${\mathit{\boldsymbol{\tau }}_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 ${\mathit{\boldsymbol{\tau }}_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}$.

4.4.   Blood flow in a representative patient-specific arterial model

In this test case we take the method for weakly imposed boundary conditions to patient-specific 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 $\Delta t = 0.005{\rm{ 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 $\rho = 1.06{\rm{ g}}{\kern 1pt} {\rm{/c}}{{\rm{m}}^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{\rm{ 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 ${\mathit{\boldsymbol{\tau }}_g}$ defined in Eq (28) is multiplied by a constant $a = 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.

5.   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 blood-artery 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.

Acknowledgments

This research was partly supported by NIH Grant No. 1R01GN135921-01. Computing resources for this work were provided by the National Center for Supercomputing Applications (NCSA) under the Blue Waters award supported by National Science Foundation (NSF). This support is gratefully acknowledged.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix A. Derivation of the boundary term

The domain-interior terms in the narrow band along the interface, $\tilde \Omega$, are converted to the boundary term on ${\Gamma _g}$ in four steps: (ⅰ) applying integration-by-parts, (ⅱ) applying the divergence theorem and assuming the boundary term to be active only on ${\Gamma _g}$, (ⅲ) neglecting the domain-interior term based on the assumption that $\Delta {\boldsymbol{\tilde v}} \approx {\bf{0}}$ in $\tilde \Omega$, and (ⅳ) substituting the expression for the fine-scales. Two modeling assumptions are employed in the second and third steps.

Assumption 1. The incremental fine-scales vanish, $\Delta {\boldsymbol{\tilde v}} = {\bf{0}}$, on the non-Dirichlet boundaries $\Gamma \backslash {\Gamma _g}$.

Assumption 2. The incremental fine-scales are non-zero only at the Dirichlet boundary ${\Gamma _g}$ and vanishingly small, $\Delta {\boldsymbol{\tilde v}} \approx {\bf{0}}$, in the domain interior.

The following example shows the application of these transformations for the viscous term.

STEP 1: Apply the integration-by parts.

STEP 2: Apply the divergence theorem to the first term.

Assuming that the incremental fine-scale field is active, $\Delta {\boldsymbol{\tilde v}} \ne {\bf{0}}$, only at the Dirichlet boundaries ${\Gamma _g}$ (Assumption 1), the boundary term on $\Gamma$ is relegated to the boundary term on ${\Gamma _g}$. Additionally, we neglect the domain-interior terms including $\Delta {\boldsymbol{\tilde v}}$ based on the assumption that the fine-scale field is vanishingly small, $\Delta {\boldsymbol{\tilde v}} \approx {\bf{0}}$, in the domain interior near the boundary $\tilde \Omega$ (Assumption 2).

STEP 3: Substituting the expression for the fine-scale in Eq (31) in the above equation, we complete the derivation of the boundary term.

Appendix B. Boundedness of the boundary term involving the inverse of the shear-rate

The boundary term in Eq (33) involving the inverse of the shear-rate ${\dot \gamma ^{{\kern 1pt} - 1}}$ can be written as

In 1D case, Eq (44) can be simplified to

where $n = - 1$ or 1.

Since the shear-rate is defined by $\dot \gamma = \sqrt 2 \left| {{\kern 1pt} \partial {\kern 1pt} v/\partial {\kern 1pt} x{\kern 1pt} } \right|$,

We realize that this term converges to zero when the shear-rate approaches to zero, $\dot \gamma \to 0$.

In 2D or 3D case, substituting $\dot \gamma = \sqrt 2 \left\| {{\kern 1pt} \mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}}){\kern 1pt} } \right\|$ in Eq (44) yields

where $\left\| \mathit{\boldsymbol{\varepsilon }} \right\|:{\kern 1pt} {\kern 1pt} = \sqrt {\mathit{\boldsymbol{\varepsilon }}:\mathit{\boldsymbol{\varepsilon }}}$. When $\dot \gamma \to 0$, ${\varepsilon _{ij}} \to 0$ and the term $\mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}})/\sqrt {\left\| {\mathit{\boldsymbol{\varepsilon }}({\boldsymbol{v}})} \right\|}$ converges to zero in the view of the L'Hopital's rule. Therefore, the boundary term in Eq (47) also converges to zero in the limit of $\dot \gamma$.