\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
Seventh Mississippi State  - UAB Conference on Differential Equations and
Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conf. 17 (2009),  pp. 171--184.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{171}
\title[\hfilneg EJDE-2009/Conf/17\hfil Finite difference methods]
{Finite difference methods for coupled flow interaction
transport models}

\author[S. McGee, P. Seshaiyer \hfil EJDE/Conf/17 \hfilneg]
{Shelly McGee, Padmanabhan Seshaiyer}  % in alphabetical order

\address{Shelly McGee\newline
Department of Mathematics, University of Findlay, Findlay,
OH 45840, USA}
\email{mcgee@findlay.edu}

\address{Padmanabhan Seshaiyer \newline
Mathematical Sciences, George Mason University, Fairfax, VA  22030, USA}
\email{pseshaiy@gmu.edu}


\thanks{Published April 15, 2009.}
\thanks{Supported by grants DMS 0813825 from the Computational
Mathematics program, NSF, \hfill\break\indent
and ARP 0212-44-C399 from the the
Advanced Research Program, Texas Higher\hfill\break\indent
 Education Coordinating Board}
\subjclass[2000]{65N30, 65N15}
\keywords{Finite difference; coupled; flow-transport}

\begin{abstract}
 Understanding chemical transport in blood flow involves coupling
 the chemical transport process with flow equations describing
 the blood and plasma in the membrane wall. In this work, we
 consider a coupled two-dimensional model with transient
 Navier-Stokes equation to model the blood flow in the vessel
 and Darcy's flow to model the plasma flow through the vessel wall.
 The advection-diffusion equation is coupled with the velocities
 from the flows in the vessel and wall, respectively to model
 the transport of the chemical. The coupled chemical transport
 equations are discretized by the finite difference method and
 the resulting system is solved using the additive Schwarz method.
 Development of the model and related analytical and numerical
 results are presented in this work.

\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{remark}{Remark}[section]

\section{Introduction}

Mathematical modelling of blood flow in the small vessels is a
challenging problem and has been of great interest to many
researchers. This problem coupled with models for the transport of
chemicals in the blood stream brings new challenges in the
solution methodology. In this paper, we present and analyze a
model that couples blood flow and chemical transport where a
discontinuous solution in the chemical concentration along one
boundary is allowed.

Specifically, we consider a mathematical model for coupling blood
flow and chemical transport in the arteries.  The Navier-Stokes
equation in two-dimensions is used to model the blood flow in an
artery and Darcy's flow is used to model the plasma flow through
the arterial wall. These respective flow equations are coupled at
the arterial wall boundary through appropriate interface
conditions.  For the chemical transport model, the
advection-diffusion equation is employed both inside the artery
and the arterial wall, and the velocities are coupled at the
arterial wall. The coupled model is discretized via the finite
difference method. In particular, we employ an explicit finite
difference method for the blood flow and an implicit finite
difference method for the advection-diffusion equation. Note that
it is also important to take into account that the concentration
across the boundary from the vessel to the wall does not have to
be continuous; we employ an iterative Schwarz type algorithm to
accomplish this.

\section{Background, models, and methods}

In this section, we will describe the models that are employed for
the blood flow, plasma flow and chemical transport on a two
dimensional domain (figure \ref{fig1}.)

\begin{figure}[htbp]\label{twoDomains}
\begin{center}
\includegraphics[width=0.3\textwidth, angle=270]{fig1} % twoDomains
\end{center}
\caption{The domain partitioned into the two domains $\Omega_f$
 and $\Omega_w$.}\label{fig1}
\end{figure}

Consider figure \ref{fig1}, where the geometry is partitioned
into two subdomains, $\Omega_f$ and $\Omega_w$, where $\Omega_f$
denotes the domain for the blood flow and $\Omega_w$ denotes the
domain for the flow of plasma in the vessel wall. The interface
between $\Omega_f$ and $\Omega_w$ is $\Gamma$, where the two
velocities are coupled as well as the chemical concentrations from
the respective advection-diffusion equatons.  Next, we describe
the equations in the respective subdomains.

\subsection{Model for blood flow}
The Navier-Stokes Equation in two-dimensions that is employed to model
the blood flow can be expressed as
\begin{gather}\label{NSmomentum}
\rho\frac{\partial \mathbf{u}}{\partial t}
+\rho(\mathbf{u}\cdot\nabla)\mathbf{u}-\nu\Delta\mathbf{u}
+\nabla p=\mathbf{f}
\\ \label{NScont}
\nabla\cdot \mathbf{u}=0.
\end{gather}
where $\mathbf{u}=\mathbf{u}(\mathbf{x},t)=(u(x,y,t),v(x,y,t))^T$
where $u(x,y,t)$ (cm/sec) is the velocity in the $x$ direction
at the point $(x,y)$ at time $t$ and $v(x,y,t)$ (cm/sec) is the
velocity in the $y$ direction at the point $(x,y)$ at time $t$,
$\rho$ is the density of the fluid (g/cm$^3$), $\nu$ is the
viscosity of the fluid (g/(cm$^2\cdot$ sec)), $p=p(x,y,t)$ is the
pressure at the point $(x,y)$ at time $t$ (g/(cm$\cdot$sec$^2$)),
and $\mathbf{f}$ is any external forces acting on the fluid at
point $(x,y)$ at time $t$ (g/(cm$^2$sec$^2$)).  Equation
(\ref{NSmomentum}) is called the \emph{momentum equation} and is
derived from the conservation of momentum which describes the
motion of the particles in the fluid. Equation (\ref{NScont}) is
called the \emph{continuity equation} which is derived from the
conservation of mass.

Note that (\ref{NSmomentum})-(\ref{NScont}) can be written in component
form as
\begin{gather}\label{NSequ}
\rho\frac{\partial u}{\partial t}+\rho\Big(\frac{\partial(u^2)}{\partial x}
+\frac{\partial (uv)}{\partial y}\Big)
-\nu\Big(\frac{\partial^2 u}{\partial x^2}
+\frac{\partial^2 u}{\partial y^2}\Big)+\frac{\partial p}{\partial x}=f_1
\\ \label{NSeqv}
\rho\frac{\partial v}{\partial t}+\rho\left(\frac{\partial
(uv)}{\partial x}+\frac{\partial (v^2)}{\partial
y}\right)-\nu\left(\frac{\partial^2 v}{\partial x^2}+
\frac{\partial^2 v}{\partial y^2}\right)+\frac{\partial
p}{\partial y}=f_2
\\\label{NScont2}
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0
\end{gather}
where $\mathbf{f}=(f_1(x,y,t),f_2(x,y,t))^T$.

We solve (\ref{NSequ})-(\ref{NScont2}) on $\Omega_f$ with boundary conditions
\begin{gather}
\mathbf{u}=u_{\rm in}(0,y,t)=U(t)(R^2-y^2)\quad\text{on }
\Omega_{f,up}\label{NSBCup}\\
\frac{\partial \mathbf{u}}{\partial x}=0\quad\text{on }
\Omega_{f,dw}\label{NSBCdw}\\
v(x,0)=0\quad\text{on }\partial\Omega_f \label{NSBCfree1}\\
\frac{\partial u(0,x)}{\partial y}=0 \quad\text{on }
\partial\Omega_f\label{NSBCfree2}\\
\mathbf{u}=0\quad\text{on }\Gamma\label{NSBCnoslip}
\end{gather}
where $u_{in}(0,y,t)$ is the prescribed inflow velocity upstream
on $\Omega_{f,up}$, $R$ is the width of the domain in the $y$ direction,
and $U(t)=(1-\cos((2\pi t+\pi)/T_p))/2$ or $U(t)$ is a constant
depending on whether the inflow is time varying or constant
\cite{griebel, karner, quarteroni}. $T_p$ is the pulse period.
 Equation (\ref{NSBCup}) is the inflow boundary condition,
and (\ref{NSBCdw}) is the outflow boundary condition downstream.
Equations (\ref{NSBCfree1}) and (\ref{NSBCfree2}) are free slip
boundary conditions. These boundary conditions allow the original
domain to be cut in half during the computations since the flow
will be symmetric about $\partial\Omega_f$. Equation
(\ref{NSBCnoslip})  is the no slip boundary condition.  An initial
condition $\mathbf{u}(\mathbf{x},0)$ is also prescribed along with
the boundary conditions.

\subsection{Model for plasma flow}

Since the blood vessel wall is porous as in saturated ground water flow, we employ Darcy's Law to describe the flow of plasma through the vessel wall \cite{karner}. This is written as
\begin{gather}\label{darcy1}
\mathbf{u}_{\rm filt}=-K\nabla p
\\ \label{darcy2}
\nabla\cdot\mathbf{u}_{\rm filt}=0
\end{gather}
where $\mathbf{u}_{\rm filt}$ is the filtration velocity of the
fluid (cm/sec), $K$ is the conductivity of the porous media
(cm$^3$sec/g), and $p$ is the pressure.  Equation (\ref{darcy2})
can be interpreted as the filtration velocity is equal over the
domain. This is a result of the porous media being saturated.
Applying the gradient operator to (\ref{darcy1}), the equation
becomes \cite{karner}
\begin{equation}
\nabla\cdot\mathbf{u}_{\rm filt}=\nabla\cdot(K\nabla p)
\end{equation}
which yields
\begin{equation}\label{laplacePressure}
0=\Delta p,
\end{equation}
if $K$ is considered to be a constant.  Equation (\ref{laplacePressure})
is solved on $\Omega_w$ with the boundary conditions given by
\begin{gather}
p=p_0(\mathbf{x}) \quad\text{on }\Gamma\label{DBCin}\\
p=p_1(\mathbf{x})\quad\text{on }\partial\Omega_{w}\label{DBCout}\\
\frac{\partial p}{\partial \mathbf{n}}=0 \quad\text{on
}\partial\Omega_{w,up}\text{ and }\partial\Omega_{w,dw}\label{DBCup}
\end{gather}
and then use the solution $p$ in (\ref{darcy1}) to solve
$\mathbf{u}_{\rm filt}$ for a known $K$.


\subsection{Chemical transport in blood flow}

The chemical transport for a flow system such as blood flow is modelled
by the advection-diffusion equation
\begin{equation}
\frac{\partial C}{\partial t}-D\Delta C+\mathbf{u}\cdot\nabla C=s
\end{equation}
with appropriate boundary conditions where $C$ is the
concentration at $(x,y)$ in $(g/cm^3)$, $D$ is the diffusivity in
$(cm^2/sec)$, $\mathbf{u}$ is the velocity of the fluid in
$(cm/sec)$, and $s$ is a  source term in $(g/(cm^3sec))$.  For the
flow structure system  described previously, two
advection-diffusion equations are needed, one to describe
concentration $C_f$ in the fluid domain, $\Omega_f$, and the other
for the concentration $C_w$ in the structure domain, $\Omega_w$.
Since the interface $\Gamma$ represents a physical barrier that is
a selectively permeable membrane (the blood vessel wall) the
concentrations $C_f$ and $C_w$ do not have to match at $\Gamma$
(see figure \ref{fig1}).  The rate at which the chemical crosses
the membrane is proportional to the gradient between the two
concentrations at the barrier $\Gamma$ \cite{quarteroni}.  The
resulting system then becomes:
\begin{equation} \label{conc_eq1}
\begin{gathered}
\frac{\partial C_f}{\partial t}-D_f\Delta C_f+\mathbf{u}_f\cdot\nabla C_f
= s_f\text{ in } \Omega_f \\
\frac{\partial C_w}{\partial t}-D_w\Delta
C_w+\mathbf{u}_w\cdot\nabla C_w
= s_w\text{ in } \Omega_w
\end{gathered}
\end{equation}
with boundary conditions
\begin{gather}
D_f\frac{\partial C_f}{\partial y}+\zeta(C_f-C_w)=0\quad\text{on }\Gamma\label{schw1}\\
D_w\frac{\partial C_w}{\partial y}+\zeta(C_w-C_f)=0\quad\text{on }\Gamma\label{schw2}\\
C_f=C_0\quad\text{on }\partial\Omega_{f,up}\\
D_f\frac{\partial C_f}{\partial x}=0 \quad\text{on }\partial\Omega_{f,dw}\\
C_w=C_0\quad\text{on }\partial\Omega_{w,up}\\
D_w\frac{\partial C_w}{\partial x}=0\quad\text{on }
\partial\Omega_{w,dw}, \label{conc_bdry}
\end{gather}
where $D_f$ is the diffusivity in $\Omega_f$, $D_w$ is the diffusivity
in $\Omega_w$,  $\zeta$ is the permeability of the membrane that $\Gamma$
represents, and is calculated from a function of the shear stress
$\tilde{\sigma}$ exerted by the blood on the wall. Initial conditions
are $C_f=C_{f0}$ and $C_w=C_{w0}$ at $t=0$.

\begin{remark} \label{rmk1} \rm
The boundary conditions in (\ref{schw1}) and (\ref{schw2}), when added
together yield
\begin{equation}
D_f\frac{\partial C_f}{\partial y}=-D_w\frac{\partial C_w}{\partial y},
\end{equation}
which says that the flux out of $\Omega_f$ is equal to the flux into
$\Omega_w$, and indeed this makes sense physically.
\end{remark}

\section{Finite Difference Methods}

Next, we describe the discretized equations for the respective
models  using the finite difference methods. First, we will
present the details of the discretization and implementation for
the flow equations.  Then we will present the discretization of
the advection-diffusion equation and details for implementing the
discontinuity in the concetrations at the arterial wall via an
Addivite Schwarz type iterative method.

\subsection{Discretization of the flow equations}

In this work, an explicit finite difference method \cite{griebel,
morton}  is employed to solve the Navier-Stokes equation.
 The discretized system for  (\ref{NSequ})-(\ref{NSeqv}) becomes
\begin{gather}\label{NSFDtimeu2}
\frac{u^{n+1}-u^n}{\Delta t}
=\frac{\nu}{\rho}\Big(\frac{\partial^2 u^n}{\partial x^2}
+\frac{\partial^2 u^n}{\partial y^2}\Big)
-\Big(\frac{\partial((u^n)^2)}{\partial x}
+\frac{\partial (uv)^n}{\partial y}\Big)
-\frac{1}{\rho}\frac{\partial p^{n+1}}{\partial x}+g^{n}_1,
\\ \label{NSFDtimev2}
\frac{v^{n+1}-v^n}{\Delta t}=\frac{\nu}{\rho}
\Big(\frac{\partial^2 v^n}{\partial x^2}
+ \frac{\partial^2 v^n}{\partial y^2}\Big)
-\Big(\frac{\partial (uv)^n}{\partial x}+\frac{\partial ((v^n)^2)}
{\partial y}\Big)-\frac{1}{\rho}\frac{\partial p^{n+1}}{\partial y}+g^n_2,
\end{gather}
where $g^n_1=\frac{1}{\rho}f_1(u^n,v^n,t_n)$ and
$g_2=\frac{1}{\rho}f_2(u^n,v^n,t_n)$.
Here $u^n=u(\mathbf{x},t_n)$, $v^n=v(\mathbf{x},t_n)$,
$t_n=t_0+n\Delta t$, $t_0$ is the initial time, $n$ is the number of
time steps from the initial time to $t_n$, and $\Delta t$ is the size
of the time step.
Solving (\ref{NSFDtimeu2}) for $u^{n+1}$ and (\ref{NSFDtimev2}) for
$v^{n+1}$ gives the equations
\begin{gather}\label{unplus1}
u^{n+1}=F^n-\Delta t\frac{1}{\rho}\frac{\partial p^{n+1}}{\partial x},
\\ \label{vnplus1}
v^{n+1}=G^n-\Delta t\frac{1}{\rho}\frac{\partial p^{n+1}}{\partial y}
\end{gather}
where
\begin{equation} \label{Fn}
\begin{aligned}
F^n&= F(u^n,v^n,t_n) \\
&= u^n+\Delta t\Big[\frac{\nu}{\rho}\Big(\frac{\partial^2 u^n}{\partial x^2}
+\frac{\partial^2 u^n}{\partial y^2}\Big)
-\Big(\frac{\partial((u^n)^2)}{\partial x}
+\frac{\partial (u^nv^n)}{\partial y}\Big)+g^n_1\Big],
\end{aligned}
\end{equation}
and
\begin{equation} \label{Gn}
\begin{aligned}
G^n&= G(u^n,v^n,t_n) \\
&= v^n+\Delta t\Big[\frac{\nu}{\rho}\Big(\frac{\partial^2 v^n}{\partial x^2}
+ \frac{\partial^2 v^n}{\partial y^2}\Big)
-\Big(\frac{\partial (u^nv^n)}{\partial x}
+\frac{\partial ((v^n)^2)}{\partial y}\Big)+g^n_2\Big].
\end{aligned}
\end{equation}
Substituting $u^{n+1}$ and $v^{n+1}$
from (\ref{unplus1})-(\ref{vnplus1}) into the discretized continuity
equation (\ref{NScont2}) we get
\begin{equation}
\frac{\partial u^{n+1}}{\partial x}+\frac{\partial v^{n+1}}{\partial y}=\frac{\partial F^n}{\partial x}-\frac{\Delta t}{\rho}\frac{\partial^2 p^{n+1}}{\partial x^2}+\frac{\partial G^n}{\partial y^2}-\frac{\Delta t}{\rho}\frac{\partial^2 p^{n+1}}{\partial y^2}=0
\end{equation}
which reduces to
\begin{equation}\label{pressurePoisson}
\frac{\partial^2 p^{n+1}}{\partial x^2}
+\frac{\partial^2 p^{n+1}}{\partial y^2}=\frac{\rho}{\Delta t}
\Big((\frac{\partial F^n}{\partial x}+\frac{\partial G^n}{\partial y}\Big).
\end{equation}
Since an initial condition is required, $u(x,y,0)=u^0$, $v(x,y,0)=v^0$,
and $p(x,y,0)=p^0$ are known, hence $F^0$ and $G^0$ can be calculated.
Thus, the right hand side of (\ref{pressurePoisson}) is known, and $p^{n+1}$
can be calculated.  Using (\ref{unplus1}) -(\ref{vnplus1}) one can
then compute $u^{n+1}$ and $v^{n+1}$.

Now the spatial derivatives are discretized also using finite
differences.   Central differences are used to approximate first
derivative terms $\frac{\partial p^{n+1}}{\partial x}$,
$\frac{\partial p^{n+1}}{\partial y}$, $\frac{\partial
F^n}{\partial x}$, and $\frac{\partial G^n}{\partial y}$.   A
staggered grid, shown in figure \ref{fig2}, is used to calculate
$u$ and $v$ because oscillations can occur in grids where $u$ and
$v$ are calculated at the same point \cite{griebel}.  Bilinear
interpolation is used to calculate $u$, $v$, and $p$ at points
other than node points \cite{press}. Notice that $\frac{\partial
p^{n+1}}{\partial x}$ is used in (\ref{unplus1}) to calculate
$u(x_i,y_{j+1/2},t_{n+1})$,  where $x_i=x_0+i\Delta x$,
$i=0,1,2,\dots ,I$, and $y_{j+1/2}=y_0+(j+1/2)\Delta y$,
$j=0,1,2,\dots ,J$, $x_0$ and $y_0$ are the smallest values in the
domain, $\Delta x$ and $\Delta y$ are the step sizes in the $x$
and $y$ direction, respectively, and $I$ and $J$ are the number of
steps in the $x$ and $y$ directions, respectively, and must be
integers. The half steps are due to the staggered grid (see figure
\ref{fig2}). Let us use the notation
$u_{i,j}^n=u(x_i,y_{j+1/2},t_n)$ and
$u_{i+1,j}^n=u(x_{i+1},y_{j+1/2},t_n)$ from now on. Simple central
differences on convection terms can cause oscillations in the
numerical solution that are not physically occurring
\cite{griebel, morton}.  The solution methodology used here
employs a donor-cell scheme.

\begin{figure}\label{staggeredGrid}
\begin{center}
\includegraphics[width=0.45\textwidth,angle=270]{fig2} % staggerdGrid.eps
\end{center}
\caption{Staggered grid for the Navier-Stokes discretization.}\label{fig2}
\end{figure}

Now the boundary conditions in (\ref{NSBCup})-(\ref{NSBCnoslip})
for the  staggered grid will be  formulated as follows. The inflow
condition given in (\ref{NSBCup}) for the velocity in the $x$
direction is
\begin{equation}\label{ICinflowU}
u_{i,j}^n=U(t_n)(R^2-y_{j+1/2}^2),\quad \text{for }i=0
\text{ and }j=1,\dots ,J.
\end{equation}
The inflow in the $y$ direction is zero so
\[
v(0,y,t)=0,
\]
but since $v$ is evaluated at the points $(x_{-1/2},y_j,t_n)$ and
$(x_{1/2},y_j,t_n)$, $v$ on the boundary is calculated by
\[
v(0,y_j,t)=v(x_{-1/2},y_j,t_n)+v(x_{1/2},y_j,t_n)=0
\]
leading to the boundary condition for the method formulated as
\[
v_{i,j}^n=-v_{i+1,j}^n\quad\text{for }i=0\text{ and }j=0,\dots ,J,
\]
using the numbering scheme shown on the grid in figure \ref{fig2}.
The outflow condition given in (\ref{NSBCdw}) is
\[
u_{i,j}^n=u_{i-1,j}^n\quad \text{for }i=I\text{ and }j=0,\dots ,J+1,
\]
and
\[
v_{i,j}^n=v_{i-1,j}^n\quad \text{for }i=I+1\text{ and }j=0,\dots ,J.
\]
The no slip and free slip boundary conditions on the remaining two
sides, $\Gamma$ and $\partial \Omega_f$, respectively,
(see figure \ref{fig1}) are
\begin{gather}
u_{i,j}^n = -u_{i,j-1}^n\quad\text{for }i=0,\dots ,I\text{ and }j=J+1,\\
v_{i,j}^n = 0\quad\text{for }i=0,\dots ,I+1\text{ and }j=J,\label{noslip2}\\
u_{i,j}^n = -u_{i,j+1}^n\quad\text{for }i=0,\dots ,I\text{ and }j=0,\\
v_{i,j}^n = 0\quad\text{for }i=0,\dots ,I+1\text{ and
}j=0\label{noslip4}.
\end{gather}


Solving the system starts with the initial condition for $u$, $v$, and $p$,
supplied by the user.  Since this is an explicit finite difference scheme,
stability conditions must be met. These conditions can be derived to be
\begin{gather*}
\frac{2\nu\Delta t}{\rho} < \Big(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}
\big)^{-1}\\
|u_{\rm max}|\Delta t < \Delta x\\
|v_{\rm max}|\Delta t < \Delta y,
\end{gather*}
where $u_{\rm max}=\max|u_{i,j}^n|$ and $v_{\rm max}=\max|v_{i,j}^n|$.
The last two conditions are the Courant-Friedrichs-Lewy (CFL)
conditions and they guarantee that no fluid particle will travel more
than $\Delta x$ or $\Delta y$ in a single time step $\Delta t$.
To ensure that all the stability conditions are met, $\Delta t$
should be chosen so that
\begin{equation}\label{stability}
\Delta t<\min\Big(\frac{\rho}{2\nu}\Big(\frac{1}{\Delta x^2}
+\frac{1}{\Delta y^2}\Big)^{-1},
\frac{\Delta x}{|u_{\rm max}|},\frac{\Delta y}{|v_{\rm max}|}\Big).
\end{equation}


After checking the stability condition and adjusting $\Delta t$ if needed,
$F^0$ and $G^0$ are calculated from $u^0$ and $v^0$. The parameter
$\gamma$ should be chosen as \cite{griebel}
\begin{equation}\label{gamma1}
\gamma=\max_{i,j}\Big(\big|\frac{u_{i,j}^n\Delta t}{\Delta x}\big|,
\big|\frac{v_{i,j}^n\Delta t}{\Delta y} \big|\Big)
\end{equation}
Using $F^0$ and $G^0$, equation (\ref{pressurePoisson}) is solved
using the iterative method, \emph{successive over-relaxation} (SOR).
The latter is used because the operations can be performed cell wise
without ever having to creating a matrix, but other methods to solve
this Poisson equation can also be used. The boundary conditions on $p$
are no flux boundary conditions except at the downstream boundary
$\Omega_{f,dn}$, which has a Dirichlet boundary condition in order
to avoid problems with uniqueness. Now $p^{n+1}$ is used in equations
(\ref{unplus1}) and (\ref{vnplus1}) to calculate $u$ and $v$ for the
next time step.  Then all the boundary conditions are updated and
the processes is repeated.

On $\Omega_w$ (see figure \ref{fig1}) the equation for the
pressure  $\Delta p=0$ or
\begin{equation}
\frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=0
\end{equation}
is solved for $p$ with Dirichlet boundary conditions on the upper
and lower boundaries and  no flux boundary conditions on the
sides.  Using a finite difference discretization via central
differences this equation can be transformed into a matrix system
that can be solved using several efficient matrix solver routines.
Then $p$ is used to solve for $u$  and $v$ with
\begin{gather}
u = -K\frac{\partial p}{\partial x},\label{uWall}\\
v = -K\frac{\partial p}{\partial y},\label{vWall}
\end{gather}
where $\frac{\partial p}{\partial x}$  and $\frac{\partial p}{\partial y}$ use central differences for first order derivatives.

\begin{remark} \label{rmk2} \rm
In the coupled system consisting of $\Omega_f$ and $\Omega_w$
shown in figure \ref{fig1}, matching flow velocities and pressure
at the  interface $\Gamma$ is required.  This is accomplished by
using the pressure at $\Gamma$ from $\Omega_f$ as the boundary
condition at $\Gamma$ on $\Omega_w$. Then after $u$ and $v$ are
calculated for $\Omega_w$, $u$ and $v$ at $\Gamma$ are used for
the boundary conditions at  $\Gamma$ in $\Omega_f$ as
\begin{gather}
u_{i,j}^n = u_{w\;i,0}^n\Delta x_f-u_{i,j-1}^n\quad\text{for }
 i=0,\dots ,I\text{ and }j=J+1,\label{noslip5}\\
v_{i,j}^n = v_{w\;i,0}^n\quad\text{for }i=0,\dots ,I+1
\text{ and}j=J,\label{noslip6}
\end{gather}
where $u_{w\,i,0}$ and $v_{w\,i,0}$ is velocity from $\Omega_w$
(which may be  interpolated values if the grids are non-matching)
and $\Delta x_f$ is the step size in the $x$ direction on $\Omega_f$.
\end{remark}

\subsection{Discretization of the coupled problem}

The time discretized system for the advection-diffusion equation
system (\ref{conc_eq1})-(\ref{conc_bdry}) can be written using an
implicit scheme in time and employing an additive Schwarz
technique \cite{gander} as:
\begin{gather}
C_f^{n+1,k+1}+\Delta t\mathcal{L}_f(C_f^{n+1,k+1})
 =\Delta t s_f^{n+1,k+1}+C_f^n\label{adeq1}\\
C_w^{n+1,k+1}+\Delta t\mathcal{L}_w(C_w^{n+1,k+1})
 =\Delta t s_w^{n+1,k+1}+C_w^n\\
D_f\frac{\partial C_f^{n+1,k+1}}{\partial y}+\zeta C_f^{n+1,k+1}
 =\zeta C_w^{n+1,k}\\
D_w\frac{\partial C_w^{n+1,k+1}}{\partial y}+\zeta C_w^{n+1,k+1}
 =\zeta C_f^{n+1,k}\\
\mathcal{B}_f(C_f^{n+1,k+1})=g_f^{n+1,k+1}\\
\mathcal{B}_w(C_w^{n+1,k+1})=g_w^{n+1,k+1}\label{adbc4}
\end{gather}
where $n$ and $n+1$ are the time steps, and $k$ and $k+1$ are the
iteration number.  $C_f^n$ and $C_w^n$ are known from the previous
time step,  $C_f^{n+1,k}$ and $C_w^{n+1,k}$ are known from the
previous iteration, and $C_f^{n+1,k+1}$ and $C_w^{n+1,k+1}$ are
the unknown values for which the system will be solved.
$\mathcal{B}_f$ and $\mathcal{B}_w$ are the boundary functions for
the boundaries of $\Omega_f$ and $\Omega_w$, respectively.  They
do not contain information about the interface $\Gamma$ shared by
both domains.

\begin{remark} \label{rmk3} \rm
The additive Schwarz method is a domain decomposition method that
allows the problem on each domain to be solved independently, then
the values common to both domains are exchanged and the solution
is calculated again.  This is repeated until the change in the
solution in successive iterations is smaller than a user supplied
tolerance \cite{gander}.  The additive Schwarz method can have a
region of overlap that in two dimensions has a positive area. In
these cases the solution on the two domains is continuous, and in
fact, the two domains could quite naturally be computed as a
single domain, but are divided to solve each domain simultaneously
on separate processors. In the case discussed here, the overlap is
only the interface $\Gamma$.
\end{remark}

Applying discretizations, the above system can be written in the matrix form
\begin{gather*}
A_1C_f^{n+1,k+1}=f_1^{n+1,k+1}+B_1C_w^{n+1,k}\\
A_2C_w^{n+1,k+1}=f_2^{n+1,k+1}+B_2C_f^{n+1,k}
\end{gather*}
where $f_1^{n+1,k+1}=s_f^{n+1,k+1}+C_f^n$, $A_1$ is the matrix
associated $\Omega_f$ and the boundaries of $\Omega_f$, $B_1$
is associated with the interface $\Gamma$ on the $\Omega_f$ side,
and the second equation is analogous for $\Omega_w$.

\begin{remark} \label{rmk4} \rm
One can check that that finite difference scheme presented for the
coupled problem is consistent and stable \cite{shelly}. One can
also derive the following error estimate: The $l_\infty$ error for
the two domain system (\ref{adeq1})-(\ref{adbc4}) can be shown to
be \cite{shelly}
\begin{equation}
E_f^{n+1}+E_w^{n+1}\leq n\Delta t \frac{\bar{K}}{1-\hat{K}}(T_f+T_w)
\end{equation}
where
\begin{gather*}
E_f^{n+1}= \max_{i,j \in \Omega_f} |{C_f}_{i,j}^{n+1}
   - C_f(x_i, y_j, t_{n+1})| \\
E_w^{n+1}= \max_{i,j \in \Omega_w} |{C_w}_{i,j}^{n+1} - C_w(x_i,
y_j, t_{n+1})|
\end{gather*}
are the respective errors between the exact solution and the finite
 difference solution at the time level $t_{n+1}$ and $T_f$ and $T_w$ are the respective truncation errors in the fluid and the wall respectively. One can show that \cite{shelly}, as $\Delta t$, $\Delta x$, and $\Delta y$ go to zero, $\bar{K}$ goes to one and $\hat{K}$ goes to zero, so that the coefficient in front of the truncation errors is going to one, and as the truncation errors go the zero the method is convergent. The details of the proof can be found in \cite{shelly}.
\end{remark}

\section{Numerical results and discussion}

In this section, we present numerical experiments that models the
methods developed in earlier sections. The method is validated by
employing exact solutions to known boundary conditions and
evaluating the magnitude of the truncation errors. Several exact
solutions are tested and the calculated finite difference solution
is compared with the exact solution.  It is shown that the error
in the method is less than the truncation error.

Consider the exact solutions
\begin{gather*}
C_f(x,y,t)=-tx(2-x)y^2\quad\text{on }0\leq x\leq 1,\ 0\leq y\leq 1\\
C_w(x,y,t)=tx(2-x)y(4-y)\quad\text{on }0\leq x\leq 1,\ 1\leq y\leq 2
\end{gather*}
that satisfy the advection-diffusion equation
\begin{equation}
\frac{\partial C}{\partial t}=D\Big(\frac{\partial^2 C}{\partial
x^2} +\frac{\partial^2 C}{\partial y^2}\Big)-u\frac{\partial
C}{\partial x}-v\frac{\partial C}{\partial y}+f
\end{equation}
where $f$ is a source term.  Here, $D$, $u$, and $v$ are all set to one
(for simplicity), and the source term on $\Omega_f$ can be calculated as
\begin{equation}
f_f=-2ty^2+2tx(2-x)-t(2-2x)y^2-2tx(2-x)y-x(2-x)y^2,
\end{equation}
and the source term on $\Omega_w$ can be calculated as
\begin{equation}
f_w=2ty(4-y)+2tx(2-x)+t(2-2x)y(4-y)+tx(2-x)(4-2y)+x(2-x)y(4-y).
\end{equation}
The boundary conditions for $C_f$ are given by
\begin{gather}
C_f(0,y,t)=0\label{eqn1}\\
\frac{\partial C_f(1,y,t)}{\partial x}=0\\
\frac{\partial C_f(x,0,t)}{\partial y}=0\\
\frac{\partial C_f(x,1,t)}{\partial y}+\zeta C_f(x,1,t)=\zeta
C_w(x,1,t)
\end{gather}
with the initial condition
\begin{equation*}
C_f(x,y,0)=0,
\end{equation*}
and the boundary conditions for $C_w$ are given by
\begin{gather*}
C_w(0,y,t)=0\\
\frac{\partial C_w(1,y,t)}{\partial x}=0\\
\frac{\partial C_w(x,2,t)}{\partial y}=0\\
\frac{\partial C_w(x,1,t)}{\partial y}+\zeta C_w(x,1,t)=\zeta
C_f(x,1,t)
\end{gather*}
with the initial condition
\begin{equation}\label{eqn2}
C_w(x,y,0)=0.
\end{equation}
To satisfy the boundary conditions for $C_f$ and $C_w$  on $y=1$,
the function $\zeta$ is chosen as $\zeta=-1/2$. One can show that
the truncation error for the advection-diffusion equation is
bounded as \cite{shelly}
\begin{equation}
|\tau|\leq \Big|C_{xxxx}\frac{\Delta x^2}{12} +C_{xxx}\frac{\Delta
x^2}{6} +C_{yyyy}\frac{\Delta y^2}{12}+C_{yyy}\frac{\Delta y^2}{6}
+C_{tt}\frac{\Delta t}{2}\Big|.
\end{equation}
For the exact solution considered, the truncation error is zero.
One therefore expects that the difference between the exact and
calculated solution will be small enough to attribute to round off
error.  After setting the code parameters to those values
specified above, the maximum absolute error between the exact
solution and the calculated solution was found to be $5.5022\times
10^{-11}$, which is indeed small enough to attribute to round off
error.

Next, we consider the exact solution
\begin{gather*}
C_f(x,y,t)=-t^2x(2-x)y^2\quad\text{on }0\leq x\leq 1,\ 0\leq y\leq 1\\
C_w(x,y,t)=t^2x(2-x)y(4-y)\quad\text{on }0\leq x\leq 1,\ 1\leq
y\leq 2,
\end{gather*}
with the boundary and initial conditions shown in (\ref{eqn1})-(\ref{eqn2}).
 The parameters $D$, $u$, and $v$ are one, and $\zeta=-1/2.$ The source
term $f_f$ in $\Omega_f$ is
\begin{equation}
f_f=-2t^2y^2+2t^2x(2-x)-t^2(2-2x)y^2-2t^2x(2-x)y-2tx(2-x)y^2
\end{equation}
and the source term in $\Omega_w$ is
\[
f_w=2t^2y(4-y)+2t^2x(2-x)+t^2(2-2x)y(4-y)+t^2x(2-x)(4-2y)+2tx(2-x)y(4-y).
\]
The maximum truncation in $\Omega_f$ is bounded by
\[
|\tau_f|\leq \frac{\Delta t}{2}\max\big|\frac{\partial^2
C_f}{\partial t^2}\big|=\Delta t\max\big|x(2-x)y^2\big|=\Delta
t=T_f
\]
and the maximum truncation in $\Omega_w$ is bounded by
\[
|\tau_w|\leq \frac{\Delta t}{2}\max
\big|\frac{\partial^2 C_w}{\partial t^2}\big|
=\Delta t\max\left|x(2-x)y(4-y)\right|=4\Delta t=T_w.
\]
So the calculated solution should be bounded by $T=T_f+T_w=5\Delta
t$. Table \ref{table} shows how the error decreased as $\Delta t$
decreased  and also shows that the error is less than the maximum
truncation error $T$.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|} \hline
$n$& $\Delta t$ & $\max|C-C_{\rm exact}|$ & $T$\\
\hline
1&1&.8224&5\\
10&.1&.0988&.5\\
100&.01&.0100&.05\\
1000&.001&.0010&.005\\
10000&.0001&.0001&.0005\\
\hline
\end{tabular}
\end{center}
\caption{Relationship between $\Delta t$,  maximum and truncation error}
 \label{table}
\end{table}
Other test cases have been considered in \cite{shelly}. The
experiments
 clearly suggested that the maximum error in the calculated solution
compared to the exact solution was less than the sum of the maximum
truncation errors on each domain as predicted theoretically.

Next, we demonstrate the performance of the numerical method for
the coupled problem. The velocities needed for the
advection-diffusion equation are obtained by solving the fluid
equations as described earlier. These equations are run until the
steady-state solution is reached. Then the $u$ and $v$ values at
the $(x,y)$ points used by the advection-diffusion equation are
calculated.  The diffusivity $D_f$ and $D_w$ is set to $1\times
10^{-5}$ as in \cite{quarteroni}.   The domain on which the
chemical transport is calculated is $3\leq x\leq10$ and $0\leq
y\leq .31$ for $\Omega_f$ and $3\leq x\leq 10$ and $.31\leq y\leq
.3414$ for $\Omega_w$.  The initial concentration is zero
everywhere except at $\Omega_{f,up}$ and $\Omega_{w,up}$, where it
is set to 100 $g/cm^3$.  The parameter $\zeta$ is calculated by
\begin{equation}
\zeta=10^{-4}(1+|\tilde{\sigma}|)
\end{equation}
where $\tilde{\sigma}$ is the stress tensor and is calculated as
in \cite{quarteroni}, with
\begin{equation}
\tilde{\sigma}=p\tilde{I}+\nu\left[\nabla \mathbf{u}+(\nabla\mathbf{u})^T\right]
\end{equation}
as in \cite{cheng} where $p$ is the pressure, $\tilde{I}$ is the
$2\times 2$ identity matrix, $\mathbf{u}$ is the velocity, and
$\nu$ is the viscosity calculated on the boundary $\Gamma$.
Figures \ref{fig3} and \ref{fig5} are the chemical concentrations
on $\Omega_f$ and $\Omega_w$, respectively, at final time.  The
pressure and velocities used in the chemical transport equations
were the solutions to the Navier-Stokes equation and Darcy's Law.
The initial concentration was set uniformly to zero and the
chemical concentration on $\partial\Omega_{f,up}$ and
$\partial\Omega_{w,up}$ was set to one hundred $g/cm^3$.  Then to
illustrate the role of pressure in allowing the chemical to cross
the arterial wall, the pressure was set to zero and the velocities
at the final time remained the same. Figures \ref{fig4} and
\ref{fig6} suggest that significantly more chemical has entered
the arterial wall in the same amount of time.

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig3} % C_fMinP_is1_8min32sec.eps 
\end{center}
\caption[Concentration on $\Omega_f$]
{Concentration on $\Omega_f$ at final time}\label{fig3}
\end{figure}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4} % C_fP_is0_8min32sec.eps
\end{center}
\caption[Concentration on $\Omega_f$ for zero pressure]
{Concentration on $\Omega_f$ at final time  for pressure uniformly zero.}
\label{fig4}
\end{figure}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig5} % C_wMinP_is1_8min32sec.eps 
\end{center}
\caption[Concentration on $\Omega_w$]
{Concentration on $\Omega_w$ at final time} \label{fig5}
\end{figure}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig6} % C_wP_is0_8min32sec.eps
\end{center}
\caption[Concentration on $\Omega_w$ for zero pressure]
{Concentration on $\Omega_w$ at final time for pressure uniformly zero.}
\label{fig6}
\end{figure}


\section{Conclusions and Future Directions}

In this paper, chemical transport was coupled with Navier-Stokes
and Darcy's Law and analyzed computationally. The solution
methodology involved an iterative algorithm that employed the
additive Schwarz method applied to two domains where the interface
had mixed boundary conditions and discontinuity in the solution at
the interface was allowed. A numerical software to simulate the
coupled process was developed that employed the finite difference
method and several numerical experiments were performed to
validate the method presented herein.

The next step in modelling the chemical transport in blood flow is
to allow the pressure in the vessel to move and change the
thickness of the wall.  We are currently studying the coupled
system involving the flow and concentration equations described in
this paper in the cylindrical coordinate system. To these we plan
to add equations that model the visco-elastic nature of the
arterial wall and study the associated fluid-structure interaction
problem. This will help us to understand the effects of the moving
wall on the diffusion of chemicals from the blood stream into the
wall.  The problem over other geometeries will also be
investigated.  Also, we have considered here that the chemical
concentration does not change the properties of the blood and
therefore does not change the velocity. We hope to consider this
effect. Finally, we plan to use the model developed in conjunction
with experimental data to estimate parameters. The current model
provides a good insight and motivation to consider all these
aspects and study the coupled fluid-structure-concentration
problem which will be the focus of our forthcoming paper.

\begin{thebibliography}{0}

\bibitem{cheng}
Cheng, C. P., Parker, D., Taylor, C. A.
\newblock ``Quantification of Wall Shear Stress in Large Blood Vessels Using
Lagrangian Interpolation Functions with Cine
Phase-Contrast Magnetic Resonance Imaging,''
\newblock \emph{Annals of Biomedical Engineering}
\newblock 30: pp 1020-1032 (2002).

\bibitem{griebel}
Griebel, M., Dornseifer, T., Neunhoeffer, T.
\newblock Numerical Simulation in Fluid Dynamics: A Practical Introduction.
\newblock Philadelphia: Society of Industrial and Applied Mathematics, 1998.

\bibitem{karner}
Karner, G., Perktold, K., Zehentner, H. P., Prosi, M.,
\newblock ``Mass transport in large arteries and through the artery wall.''
\newblock In \emph{Intra and Extracorporeal Cardiovascular Fluid Dynamics, Volume 2-Fluid Structure Interactions},
\newblock ed. P. Verdonck and K. Pertold,
\newblock 209-247.
\newblock Southampton, Boston: WIT Press, 2000.

\bibitem{shelly}
McGee, S. M.,
\newblock \emph{Computational modeling of chemical transport in flow structure interactions.}
\newblock Doctoral disseration,
\newblock Texas Tech University, Lubbock, 2007.

\bibitem{morton}
Morton, K. W., Mayers, D. F.
\newblock \textit{Numerical Solutions of Partial Differential Equations.}
\newblock Cambridge: Cambridge University Press, 1994.

\bibitem{press}
Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.
\newblock \textit{Numerical Recipes in C++}.
\newblock Cambridge: Cambridge University Press, 2002.

\bibitem{quarteroni}
Quateroni, A., Veneziani, A., Zunino, P.
\newblock ``A domain decomposition method for advection-diffusion processes with application to blood solutes,''
\newblock \emph{SIAM Journal of Scientific Computing}
\newblock 23: pp 1959-1980 (2002).

\bibitem{gander}
St-Cyr, A., Gander, M. J., Thomas, S. J.
\newblock ``Optimized Restricted Additive Schwarz Methods.''
\newblock In \emph{16th International Conference on Domain Decomposition.,}
\newblock New York: On-line, 2004.

\end{thebibliography}
\end{document}