$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

For arterial blood flow, backflow stabilization has been proposed previously to deal with truncated domains where the entire artery is not modeled explicitly [17, 31]; for these types of problems, letting $\mathbf{t}=\mathbf{0}$ or prescribing a constant pressure at the outflow boundary may not prevent flow reversals that compromise convergence of an analysis. Instead, these authors proposed a velocity-dependent traction boundary condition, $\mathbf{t}=\beta\rho\left(\mathbf{v}\otimes\mathbf{v}\right)\cdot\mathbf{n}$ with a tensile normal component, that counters the backflow (only when $v_{n}<0$ ). Here, $\beta$ is a non-dimensional user-defined parameter; a value of $\beta=0$ turns off this feature, while a value of $\beta=1$ generally shows good numerical performance. We adapt this previously proposed formulation by letting the normal component of the viscous traction be given by $$\begin{equation} t_{n}^{\tau}=\begin{cases} \beta\rho_{r}v_{n}^{2} & v_{n}<0\\ 0 & v_{n}\ge0 \end{cases}\,.\label{eq:BFS-normal-traction} \end{equation}$$ The choice of $\rho_{r}$ in lieu of $\rho$ is for convenience, to avoid the dependence of $\rho$ on $J$ (which is negligible for nearly incompressible flow). Then, the contribution of this traction to the virtual external work $\delta W_{ext}$ is $$\begin{equation} \delta G=\int_{\partial\Omega}\delta\mathbf{v}\cdot t_{n}^{\tau}\mathbf{n}\,da\,.\label{eq:BFSckflow-virtual-work} \end{equation}$$ The linearization of $\delta G$ along an increment $\Delta\mathbf{v}$ in the velocity is given by $$\begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\int_{\partial\Omega}\delta\mathbf{v}\cdot\mathbf{K}_{n}\cdot\Delta\mathbf{v}\,da\,,\label{eq:BFS-linearization} \end{equation}$$ where $$\begin{equation} \mathbf{K}_{n}=\begin{cases} 2\beta\rho_{r}v_{n}\left(\mathbf{n}\otimes\mathbf{n}\right) & v_{n}<0\\ \mathbf{0} & v_{n}\ge0 \end{cases}\,.\label{eq:BFS-stiffness} \end{equation}$$ The discretized form of $\delta G$ is $$\begin{equation} \delta G=\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}^{n}\,,\quad\mathbf{f}_{a}^{n}=\int_{\partial\Omega}N_{a}t_{n}^{\tau}\mathbf{n}\,da\,,\label{eq:BFS-discretized-work} \end{equation}$$ whereas the discretized form of $D\delta G\left[\Delta\mathbf{v}\right]$ is $$\begin{equation} D\delta G\left[\Delta\mathbf{v}\right]=\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{n}\cdot\Delta\mathbf{v}_{b}\,,\quad\mathbf{K}_{ab}^{n}=\int_{\partial\Omega}N_{a}N_{b}\mathbf{K}_{n}\,da\,.\label{eq:BFS-discretized-tangent} \end{equation}$$