# A new pseudo-spectral methodology without numerical diffusion for conducting dye simulations and particle residence time calculations

Dye experimentation is a widely used method in experimental fluid mechanics for flow analysis or for the study of the transport of particles within a fluid. This technique is particularly useful in biomedical diagnostic applications ranging from hemodynamic analysis of cardiovascular systems to ocular circulation. However, simulating dyes governed by convection-diffusion partial differential equations (PDEs) can also be a useful post-processing analysis approach for computational fluid dynamics (CFD) applications. Such simulations can be used to identify the relative significance of different spatial subregions in particular time intervals of interest in an unsteady flow field. Additionally, dye evolution is closely related to non-discrete particle residence time (PRT) calculations that are governed by similar PDEs. PRT is a widely used metric for various fluid dynamics applications (e.g., environmental fluids, biological flows) and is a well-accepted biomarker for cardiovascular diseases since it is linked to thrombus formation. This contribution introduces a pseudo-spectral method based on Fourier continuation (FC) for conducting dye simulations and non-discrete particle residence time calculations without numerical diffusion errors. Convergence and error analyses are performed with both manufactured and analytical solutions. The methodology is applied to three distinct physical/physiological cases: 1) flow over a two-dimensional (2D) cavity; 2) pulsatile flow in a simplified partially-grafted aortic dissection model; and 3) non-Newtonian blood flow in a Fontan graft. Although velocity data is provided in this work by numerical simulation, the proposed approach can also be applied to velocity data collected through experimental techniques such as from particle image velocimetry.

## Authors

• 2 publications
• 3 publications
• 2 publications
01/28/2021

### Neural Particle Image Velocimetry

In the past decades, great progress has been made in the field of optica...
01/21/2021

### On the Use of Computational Fluid Dynamics (CFD) Modelling to Design Improved Dry Powder Inhalers

Purpose: Computational Fluid Dynamics (CFD) simulations are performed to...
11/01/2017

### A Coupled Lattice Boltzmann Method and Discrete Element Method for Discrete Particle Simulations of Particulate Flows

Discrete particle simulations are widely used to study large-scale parti...
07/15/2019

### Hydrodynamic Simulations using GPGPU Architectures

Simulating the flow of different fluids can be a highly computational in...
05/27/2020

### A numerical algorithm for Fuchsian equations and fluid flows on cosmological spacetimes

We consider a class of Fuchsian equations that, for instance, describes ...
03/11/2021

### Frame-independent vector-cloud neural network for nonlocal constitutive modelling on arbitrary grids

Constitutive models are widely used for modelling complex systems in sci...
02/03/2020

### Diffusion bridges for stochastic Hamiltonian systems with applications to shape analysis

Stochastically evolving geometric systems are studied in geometric mecha...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The tracing of dye injections constitutes a fluid visualization technique for the tracking of flows. Dye experimentation is widely used in experimental fluid mechanics for flow analysis or for the study of the transport of particles within a fluid bradshaw2016experimental. Such a method is particularly useful in biomedical diagnostic applications ranging from blood flow analysis in the cardiovascular system to ocular circulation. Although dye tracing is most commonly known as an experimental approach, it can also serve as a useful post-processing analysis tool in computational fluid dynamics (CFD) studies kim2004simulated. Numerical dye studies can facilitate the identification of the relative significance of spatial subregions within particular time intervals of an unsteady flow field. For example, dye simulations can be applied to the aorta in order to visualize the time histories of the effects of the pulsatile velocity fields on injected medication as well as to investigate the effects of previous cardiac cycles on the distribution of a drug thromb5; zhang2019lattice; kim2004simulated. Mathematically speaking, dye simulations are also related to particle residence time (PRT) calculations, both of which can be formulated as convection-diffusion problems patankar2018numerical; jozsa; kim2004simulated. PRT is an important metric for a variety of applications including those found in environmental fluid dynamics maxwell2019exploring; shen2004calculating, biological flows moran1992short, and hemodynamics reza2019critical; marsden1; marsden2. In particular, PRT is now a well-accepted biomarker for cardiovascular diseases (e.g., stroke, myocardial infarction and embolism) since it is inextricably linked to thrombus (blood clot) formation thromb1; thromb2; thromb3; thromb4; thromb5.

While significant effort has been made for the development of various PRT methods (e.g., discrete/non-discrete formulations jozsa; kim2004simulated; marsden1; marsden2, Eulerian/Lagrangian formulations maxwell2019exploring; thromb5; reza2019critical), little effort has been put towards conducting dye simulations with non-zero diffusion coefficients. Furthermore, to the best of the authors’ knowledge, there is currently no numerical methodology in use for dye simulations that incur little to no numerical diffusion errors (sometimes known as “pollution” errors). This manuscript introduces a new pseudo-spectral method based on a Fourier continuation (FC) approach for solving the two- and three-dimensional (2D, 3D) convection-diffusion partial differential equation (PDE) systems that govern dye movement and particle residence time evolution without such numerical diffusion errors.

FC-based methods broaden the applicability of Fourier-based spectral methods by enabling Fourier series representations of non-periodic functions while avoiding the well-known “Gibb’s phenomenon.” PDE solvers based on such an approach maintain the desirable qualities of spectrally-accurate numerical solvers for time-dependent systems: they can produce accurate solutions by means of relatively coarse discretizations and, importantly, they carry minimal numerical diffusion or dispersion errors (manifesting as artificial amplitude decay or period elongation). FC-based PDE solvers have been successfully constructed for a variety of physical equations including classical wave equations lyonbrunoII; brunoprieto, non-linear Burgers systems bruno_jimenez, Euler equations shz; amlanibhat, Navier-Stokes equations albinbruno; cubillos; fontanabruno, radiative transfer equations gaggiolibruno, Navier elastodynamics equations amlanibruno; amlanicarlos and 1D fluid-structure hemodynamics equations amlanipahlevan

. The numerical methodology introduced in this work represents a fully-realized 3D convection-diffusion solver based on this Fast Fourier Transform (FFT)-speed Fourier continuation approximation procedure for both Dirichlet and Neumann boundary conditions.

Convergence and error analyses are presented by considering both the method of manufactured solutions as well as analytical solutions (i.e., based on Gaussians). The performance of the proposed FC methodology is also demonstrated through three physical/physiological applications: 1) flow over a 2D cavity; 2) pulsatile flow in a simplified partially-grafted aortic dissection model; and 3) non-Newtonian blood flow in a Fontan graft.

## 2 A high-order method for dye simulations and non-discrete residence time calculations

This section introduces a new high-order, pseudo-spectral numerical method for dye simulations and non-discrete PRT calculations governed by forced and unforced convection-diffusion equations with Dirichlet or Neumann boundary conditions. The methodology produces solutions based on an FFT-speed Fourier continuation approach lyonbrunoI; albinbruno; amlanibruno; amlanipahlevan

for the accurate trigonometric interpolation of non-periodic functions, ultimately enabling high-order spatial accuracy with coarse discretizations using mild Courant-Friedrichs-Lewy (CFL) constraints on time integration. Most importantly, such a (spectral) methodology has been demonstrated to faithfully preserve the dispersion/diffusion characteristics of the underlying continuous problems (errors do not compound over distance and time, which may not be true for finite difference or finite element-based methods

kalinowska2007; sankar1998; mullen; amlanipahlevan; amlanibruno). These properties are particularly favorable for the accurate assessment of the evolution of an injected dye or drug that may circulate over long times in physiological blood flow configurations (whose corresponding velocity fields, used as input here, are provided by any appropriate fluid-structure solver, e.g., the immersed boundary-lattice Boltzmann method described later.

### 2.1 Governing equations

Given an (incompressible) fluid velocity field on a domain with fluid boundary , the corresponding evolution of the concentration of a dye kim2004simulated; marsden1 or of a non-discrete blood flow particle residence time formulation jozsa; reza2019critical; marsden2 can be generally governed by the convection-diffusion system

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂c∂t(x,t)+v(x,t)⋅∇c(x,t)−∇⋅D(x,t)∇c(x,t)=h(x,t),x∈Ω, t≥0,∇c(x,t)⋅n=hN(x,t),x∈∂ΩN, t≥0,c(x,t)=hD(x,t),x∈∂ΩD, t≥0, (1)

where ; is the unknown dye concentration (or residence time, as explained in what follows); is the diffusion coefficient; and are generic forcing functions defined, respectively, on the domain , on the Neumann boundaries ( is the outward unit normal) and on the Dirichlet boundaries . This most general formulation of a convection-diffusion system is treated by the FC solver introduced in the next section; the dye and particle residence time models of interest are merely particular choices of boundaries and forcing functions (described in what follows).

#### 2.1.1 Dye simulation problem formulation

Defining to be the flow inlet boundary (where a Dirichlet condition is imposed) and all other boundaries (where a Neumann condition is imposed), the evolution of the concentration of an injected dye (or drug) along a subset of the inlet boundary is governed by Equation (1) with kim2004simulated; marsden1

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂c∂t(x,t)+v(x,t)⋅∇c(x,t)−∇⋅D(x,t)∇c(x,t)=0,x∈Ω, t≥0,∇c(x,t)⋅n=0,x∈∂ΩN,c(x,t)=1,x∈Γdye⊂∂ΩD,c(x,t)=0,x∈∂ΩD∖Γdye, (2)

where the Dirichlet value of unity on indicates points at which a concentration is injected (e.g., at a tip of a catheter rather than the entire inlet).

#### 2.1.2 Non-discrete particle residence time problem formulation

The convection-diffusion system of Equation (1) also enables a “non-discrete” (post-processing) formulation for PRT calculations jozsa; reza2019critical; marsden2 of a given flow velocity field (which can be provided either by numerical simulation or experimentation such as through particle image velocimetry willert). Again, defining to be the flow inlet boundary and all other boundaries, the particle residence time is governed by Equation (1) with jozsa

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩∂τ∂t(x,t)+v(x,t)⋅∇τ(x,t)−∇⋅D(x,t)∇τ(x,t)=1,x∈Ω, t≥0,∇τ(x,t)⋅n=0,x∈∂ΩNc(x,t)=0,x∈∂ΩD, (3)

where, for notational purposes, for PRT calculations (usually given in units of seconds). The Dirichlet value of zero at the inlet follows from the particular configurations of interest (i.e., hemodynamics applications): since fluid flows downstream from the inlet, particles along the inlet have not spent time inside the domain For some PRT solvers used for hemodynamics applications marsden1; marsden2, the diffusion is often taken to be zero (no such assumption is required for the proposed methodology of this paper). For those that have considered mass diffusivity kim2004simulated; reza2019critical, no analysis as to numerical diffusion errors have been presented (finite difference or finite element methods are well-known to suffer from numerical diffusion kalinowska2007; sankar1998; mullen; amlanipahlevan; amlanibruno, including for convection-diffusion equations sankar1998).

### 2.2 A pseudo-spectral Fourier continuation methodology without numerical diffusion

#### 2.2.1 Spatial discretization

For the numerical treatment of the spatial variables and directional derivatives of Equation (1), one can consider discrete point values of a given smooth function (defined on the unit interval without loss of generality) through a uniform discretization of size , i.e., . The FC method constructs a fast-convergent Fourier series (or trigonometric interpolant) (where is only slightly larger than ) that is given by

 ccont=W∑k=−Wake2πikxbs.t.ccont(xi)=c(xi), i=0,...,N−1, (4)

where is the corresponding bandwidth for an -sized truncated Fourier series (here, is defined as the number of discrete points added in the interval , i.e., ). The continued (or “extended”) function renders the original periodic, discretely approximating to very high-order in but demonstrating periodicity on the slightly larger . Directional spatial derivatives on a dimension of in Equation (1) can then be calculated by exact termwise differentiation of the trigonometric series given by Equation (4) as

 ∂ccont∂x(x)=W∑k=−W(2πikb)ake2πikxb. (5)

By restricting the domain of to the original unit interval, the numerical derivatives of are hence approximated to high-order. The overall errors are found only in the production of Equation (4), from where the calculation of the derivative Equation(5) is easily facilitated by a Fast Fourier Transform (FFT). This overall idea can also be easily extended to any general interval via affine transformations lyonbrunoI.

In a most simple treatment brunohan; boyd, the coefficients of Equation (4) can be found through the solution of a least squares system given by

 minakN−1∑i=0∣∣ccont(xi)−c(xi)∣∣2, (6)

which can be solved through a Singular Value Decomposition (SVD). However, for a PDE that evolves in time such as the convection-diffusion equations for dye simulations, time integration can incur significant computational expense (where the SVDs must be applied at each timestep in multiple spatial dimensions). For tractability, this work employs an accelerated continuation technique known as FC(Gram)

lyonbrunoI; lyonbrunoII; albinbruno. Relying on the use of very small (e.g., 5) numbers of function values near the left () and right (

) endpoints, the FC(Gram) algorithm projects such values onto an interpolating (Gram) polynomial basis. Fourier continuations are precomputed for each (significantly reduced) basis vector through solving the corresponding least squares formulation similar to Equation (

6) by a high-precision or symbolic SVD. That is, via a subset of given function values on small numbers and of matching points and at the leftmost and rightmost sub-intervals of , one can produce the discrete periodic extension of size by projecting these end values onto a Gram polynomial basis up to degree (or ). This constructs an interpolating polynomial whose Fourier continuations can be precomputed, where orthogonality is enforced by the scalar product naturally defined by the discretization points. This ultimately forms a (precomputed) “basis” of continued functions that can represent any function and provide a smooth transition from to over the slightly larger interval .

The solver introduced in this work invokes a particular “blend-to-zero” variation of FC(Gram) albinbruno; amlanibruno; amlanipahlevan. Instead of continuing to the opposite endpoint of the function, the blending algorithm transitions the matching points on the left and the right smoothly to values of zero at the discrete points and . The sum of the two leftward and rightward continuations constructs the complete (discretely-periodic) Fourier continuation. Each blending transition is precomputed on each element of the Gram polynomial basis (that interpolates the handful of endpoint values of ) by solving the corresponding SVD minimization. At the left, for example, this Gram basis can be constructed by orthogonalizing a Vandermonde matrix

 V=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝1x0(x0)2...(x0)dℓ−11x1(x1)2...(x1)dℓ−1⋮⋮⋮⋮⋮1xd−1(xd−1)2...(xdℓ−1)dℓ−1⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (7)

of point values at the discrete points of the monomials (a substitution of gives the equivalent representation for the rightward matching problem). An orthogonalization can then be performed by employing Gram-Schmidt () in high-precision (see albinbruno; amlanibruno for details). This ultimately produces projection operators and (on the left and the right, respectively) in order to determine the coefficients in the polynomial bases for interpolating the matching endpoints of the function . Defining vectors of such sets of points for the left and right by

 (8)

the overall FC operation can hence be expressed as

 (9)

where contains the discrete point values of ; is a vector of the discretely-periodic function values;

is the identity matrix; and

contain the corresponding values that blend, to zero, the left and the right Gram basis vectors. The columns of , (containing the point values of each element of the corresponding Gram basis) differ only by column-ordering for (and hence by row). Further technical details on the blend-to-zero FC(Gram) algorithm have been provided previously lyonbrunoI; albinbruno; amlanibruno.

For the generalized system given by Equation (1), Neumann (derivative) conditions are often invoked for residence time and dye simulations at boundaries away from the inlet. The methodology outlined above requires function values at the absolute left and right physical endpoints (as part of the matching sets of points), which is naturally satisfied by a Dirichlet boundary-value problem. If a Neumann boundary value is specified (for example, at the right endpoint as ), it is necessary to match a given derivative value in addition to given function values contained in the rest of the matching interval (i.e., ). This can be accomplished by following a similar procedure to the classic FC(Gram) outlined above, as described previously amlanibruno: a modified polynomial interpolant can be introduced to match the derivative value at the endpoint as well as the function values at . Such an interpolant can be obtained by orthonormalizing, instead of the columns of Equation (7), the corresponding columns of a modified Vandermonde matrix given by

 Vmod=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝012x0...(dℓ−1)(x0)dℓ−21x1(x1)2...(x1)dℓ−1⋮⋮⋮⋮⋮1xdℓ−1(xdℓ−1)2...(xdℓ−1)dℓ−1⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (10)

by means of a QR-decomposition that similarly produces the operators

 Vmod=QmodRmod. (11)

In order to use the same precomputed Fourier continuation basis data produced for the Dirichlet case, a reconstruction of the coefficients from the Neumann case to the original Gram polynomial basis can be simply obtained by replacing the operator (resp. ) in Equation (9) with a modified operator (resp. ) that directly performs such a (re-)projection. The new operator can be found by first solving for new coefficients via the expression

 Vmoda=(c0,c1,...,cd−2,c′(xd−1))T (12)

which, from the QR-decomposition in Equation (11), can be solved as

 a=R−1modQTmod(c0,c1,...,cd−2,c′(xd−1))T. (13)

From the original decomposition of the original Vandermonde matrix of Equation (7), the corresponding coefficients are similarly given by

 Va=QRa(c0,c1,...,cd−2,cd−1)T. (14)

Substitution of Equation (13) into Equation (14) yields the expression

 QRR−1modQTmod(c0,c1,...,cd−2,c′(xd−1))T=(c0,c1,...,cd−2,cd−1)T.

Defining further gives

 ¯QT(c0,c1,...,cd−2,c′(xd−1))T=QT(c0,c1,...,cd−2,cd−1)T. (15)

Hence the new block matrix form of the continuation with the modified operators are by

 ccont=[cAℓ¯QTℓ¯cℓ+ArQTrcr]=[cPℓ¯cℓ+Prcr,] (16)

where and where are the same original blend-to-zero operators from the Dirichlet case. Such a formulation enables mixing and matching of left-right boundary conditions (e.g., Neumann-Dirichlet, Neumann-Neumann), where a similar procedure can be performed to find the corresponding projection operator for Neumann boundary conditions at the right endpoint.

In summary, the FC(Gram) algorithm appends values to a given discretized function in order to form a periodic extension in that transitions smoothly from back to (or, in the case of known Neumann boundaries, back to, e.g., ). An illustrative example is provided in Figure 1, where a non-periodic function is originally discretized on and then translated by a distance with the subsequent interval filled-in by the sum of leftward and rightward blend-to-zero continuations. The resulting continued vector can be interpreted as a set of discrete values of a smooth and periodic function that can be approximated to high-order via FFT on an interval of size (and hence accurate termwise derivative representations for Equation (1)). As has been suggested previously amlanibruno; amlanipahlevan, a small number of matching points with a periodic extension comprised of points are used for all simulations presented hereafter. The matrices and are computed only once (before numerically solving the governing system) and stored in file for use each time the FC procedure is invoked.

#### 2.2.2 Temporal integration

An explicit Adams-Bashforth scheme of order four (AB4) is employed for marching forward-in-time, similarly to other FC-based solvers albinbruno; amlanibruno; amlanipahlevan; amlanibhat. Other explicit timesteppers that provide adequate regions of stability bashforth; rk4musa can also be employed, although timesteppers such as the fourth-order Runge-Kutta method may entail four evaluations of the right-hand-side that is given (from Equation (1)) by

 ∂c∂t(x,t)=−v(x,t)⋅∇c(x,t)+∇⋅D(x,t)∇c(x,t)+h(x,t)=F(x,t). (17)

Multiple evaluations of at each timestep for RK4 can become rather costly for 2D and 3D simulations (where the spatial derivatives must be computed at each evaluation), and hence an AB4 method is used instead. For a uniform timestep and a given time , the corresponding equation for integrating to is given by

 ∂c∂t(x,t+Δt)=c(x,t)+Δt24(55F(x,t)−59F(x,t−Δt)+37F(x,t−2Δt)−9F(x,t−3Δt)). (18)

Dirichlet boundary conditions can then imposed after each timestep by direct injection, and Neumann boundary conditions are implemented naturally when the modified Fourier continuation projection operator is invoked to compute a spatial derivative at . For all simulations in this paper, a timestep of

 Δt≤√221max∥v(x,t)∥minΔxi (19)

has provided absolute stability in all cases, including runs of millions of timesteps.

Similarly to other spectral solvers, a spatial filter is employed at each timestep above to control the growth of errors in unresolved modes albinbruno; amlanibruno; amlanipahlevan. Applied in Fourier space on a function with Fourier coefficients , such a filter is given by

 N2∑k=−N2^gkexp(ikx)⟶N2∑k=−N2exp(−α(2k/N)p)^gkexp(ikx), (20)

where is a positive integer that determines the rate of decay of the coefficients, and where the real valued determines the level of suppression such that the highest-frequency modes are multiplied by For all numerical in this paper, values of are employed similarly to other FC-based solvers (further discussion on filtering can be found elsewhere amlanibruno; ellingthesis; amlanipahlevan).

### 2.3 Fluid-solid interaction solver: an immersed boundary-lattice Boltzmann method

The input fluid velocity field of Equation (1) can be provided by any suitable flow solver. For the hemodynamics applications of this work, such velocity vectors are produced by a non-Newtonian fluid-structure solver based on a lattice Boltzmann method for the fluid (which has been proven particularly useful for hemodynamics applications randles; ames2020multi; wei2020; wu2017lattice). This method employs simplified kinetic equations in conjunction with a modified molecular-dynamics approach, where the motions of particles on a regular lattice are enforced by distribution functions that provide mass and momentum conservation. For a position and a time , the distribution functions of this work are governed by the Bhatnagar-Gross-Krook collision model bhatnagar through the expressions

 zfi(x+eiΔt,t+Δt)−fi(x,t)=−1T[fi(x,t)−feqi(x,t)]+ΔtFi, (21)

where are distribution functions of the particles in phase space ( is the corresponding equilibrium distribution); are discrete velocities; is the temporal discretization; and is a (non-dimensionalized) relaxation time that is related to fluid viscosity. For a spatial discretization , the corresponding pressure and flow velocity are given respectively by

 P=v2s∑ifi, (22)

and

 v=1ρ∑ieifi+12fΔt, (23)

where is the lattice sound speed. In all the case studies presented herein, a so-called D2Q9 stencil is used, yielding discrete velocity vectors. Further details on the lattice Boltzmann method (LBM) can be found in literature mohamad2011lattice; chen1998lattice.

For problem formulations including compliant (elastic) walls (e.g., the simplified partially-grafted aortic dissection model described later), a corresponding structure equation can be considered in Lagrangian coordinates as

 ρsh∂2X∂t2=∂∂s[Eh(1−(∂X∂s⋅∂X∂s)−1/2)∂X∂s−∂∂s(EI∂2X∂s2)]+FL, (24)

where is the density of the solid wall; is the thickness; is the Lagrangian coordinate along the solid wall; is the position (deformation) of the solid wall; is the Lagrangian force exerted on the solid wall by the fluid; and and are the stretching stiffness and bending stiffness, respectively. The boundary condition of the solid wall at a free end is given by

 1−(∂X∂s⋅∂X∂s)−1/2=0,   ∂2X∂s2=(0,0),   ∂3X∂s3=(0,0), (25)

and at a fixed end (simply supported) is given by

 X=Xo,∂X∂s=(−1,0). (26)

Solving the solid equations via a finite element method (FEM), the immersed boundary (IB) method goldstein1993modeling; ames2020multi is used to ultimately couple the solid solution with the LBM fluid solver huang2007simulation; hua2014dynamics. This method has been extensively used to simulate fluid-structure interaction (FSI) problems in cardiovascular biomechanics  peskin2002immersed; mittal2005immersed. The reference quantities used to non-dimensionalize both the fluid and solid formulations are given by the density ; velocity (where is the inlet area); and length ( is the inlet diameter). The subsequent non-dimensional physical parameters are given by the Reynolds number ; the Womersley number ; the bending coefficient ; the tension coefficient ; and the mass ratio of the solid wall to the fluid given by .

## 3 Convergence & error analysis

### 3.1 Method of manufactured solutions for verifying implementation, accuracy & convergence

For the proposed FC-based methodology, the method of manufactured solutions amlanibruno; amlanipahlevan; roache; roy; vedovoto; mmsraghu; mmsraghu2 (MMS) has been employed in order to verify both its numerical accuracy as well as its code implementation. Such a method postulates a (sufficiently smooth) solution of the convection-diffusion PDE and, upon substitution, incorporates the corresponding right-hand side and boundary conditions as non-trivial forcing terms. The resulting system enables the proposed function to be an exact solution of a forced convection-diffusion system. For example, one can postulate a solution to as

 c(x,t)=cos(2πf1(x1−κ1t))cos(2πf2(x2−κ2t)), (27)

which, upon substitution into the PDE of Equation (1), yields a corresponding right-hand-side as

 h(x,t) = 2π[κ1f1cos(2πf2(x2−κ2t))sin(2πf1(x1−κ1t)) + κ2f2cos(2πf1(x1−κ1t))sin(2πf2(x2−κ2t))] − 2πv(x,t)⋅(f1cos(2πf2(x2−κ2t))sin(2πf1(x1−κ1t))f2cos(2πf1(x1−κ1t))sin(2πf2(x2−κ2t))).

Assuming a manufactured closed-form velocity field given by

 v=(cos(2πx)cos(12πt),sin(2πy)cos(12πt))T,

and the parameters Hz, m/s, m/s and m/s, Figure 2 presents a snapshot of the corresponding numerical solution given by (27

) at an arbitrarily chosen moment in time in a domain defined by

.

Applying the FC solver on a series of discretizations corresponding to integer divisions of (where all timesteps are chosen small enough so that errors are dominated by spatial discretizations), Figure 3 presents the maximum absolute errors, over all space and all time (up to timesteps), for both Dirichlet boundary conditions (Figure 3, left) and Neumann boundary conditions (Figure 3, right). The overlaid slopes in these plots illustrate the expected fifth-order accuracy of the operators from the FC parameters employed in this work.

### 3.2 Assessment of numerical diffusion errors: comparison with natural analytical solutions

Numerical diffusion (or “pollution”) errors of the proposed FC methodology can be studied by considering test problems for (dye) concentrations that travel over arbitrarily long distances. For a 2D domain of size m m with a constant velocity field and diffusion constant , the transport of a Gaussian pulse (ball) of unit height centered at is given by

 c(x,t)=aa+4texp[−(x1−0.1−w1t)2D(a+4t)−(x2−0.1−w2t)2D(a+4t)] (28)

and, for zero diffusion (), the solution is given by

 c(x,t)=exp[−(x1−0.1−w1t)2b1−(x2−0.1−w2t)2b2], (29)

where and are arbitrary constants.

Together with their corresponding initial conditions at , it is straightforward to verify that Equations (28) and (29) are analytical (i.e., not manufactured) solutions to the unforced convection-diffusion equation for zero and constant diffusion (respectively) of Equation (1), i.e.,

 ∂c∂t(x,t)+v0⋅∇c(x,t)−∇⋅D∇c(x,t)=0.

Snapshots of the corresponding solutions for , m/s, m/s and are presented in Figure 4 for zero physical diffusion (left; multimedia view) and non-zero physical diffusion ( m/s, right; multimedia view). Note that the constants and for the zero-diffusion case are in terms of the diffusive case coefficients in order to affect the same initial Gaussian pulse width for comparative purposes. Figure 5 presents the corresponding temporal evolution of the peak dye concentration () as well as total dye volume () produced by the FC-based solver using a spatial discretization. The corresponding analytical results marked in the figures indicate excellent agreement between analytical and numerical simulations for both the zero physical diffusion and non-zero physical diffusion cases (a weak indication of minimal numerical diffusion). The dotted lines in Figure 5 (right) demarcate the time at which the dye solution begins to exit the domain, which is expectedly sooner for the physically-diffusive case whose dye concentration has spread wider than in the non-diffusive case.

Simulation of dye propagation over long distances can be achieved by increasing the distance traveled while maintaining a fixed resolution over a characteristic length . This is analogous to wave-based solutions, for which the equivalent problem formulation is achieved by increasing the number of wavelengths across a domain while maintaining a fixed number of points-per-wavelength. Figure 6 presents a one-dimensional illustration of the problem formulation for studying numerical diffusion errors. Defining to be the effective initial width (diameter) of the Gaussian ball solution (Equation (29)), and defining to be the number (density) of points discretizing the length , an increase in corresponds to a systematic increase in the distance traveled by the peak dye concentration while maintaining this fixed resolution (i.e., a distance discretized by ).

With a timestep chosen small enough so that numerical errors are dominated by those arising from the spatial discretization, Figure 7 presents maximum numerical errors over all time and space (as functions of the problem size ) of the FC-based simulations for both the zero-diffusion and physical-diffusion cases. Each Gaussian ball is discretized by a fixed density of or points per effective diameter (). The figures are overlaid with the expected behavior of certain numerically diffusive methods kalinowska2007; amlanipahlevan; babuska; mullen such as finite differences (including for convection-diffusion equations sankar1998) or finite elements where, even for small problem sizes , larger and larger densities would be necessary in order to produce reasonable accuracies. By contrast, the numerical accuracy resulting from the FC-based solutions (for both the physically-diffusive and non-diffusive cases) remains essentially constant with increasing . The errors are effectively independent of Gaussian ball width (and distance traveled) when the density of spatial points discretizing the width is kept constant—a strong indication that the proposed FC-based methodology incurs no numerical diffusion errors.

## 4 Physical case studies

This section presents three different physical case studies: 1) a classical 2D example of steady flow over a cavity marsden2; 2) a 3D-axisymmetric example of pulsatile blood flow in a simplified partially-grafted aortic dissection (chosen for its combination of rigid and elastic walls in solid-fluid interactions); and 3) a 3D example of a non-Newtonian fluid in a Fontan graft wei2020significance .

### 4.1 Flow over a cavity (2D)

The 2D computational domain of a cavity underneath a channel cross-flow is given in Figure 8 (left) and corresponds to an example previously proposed for residence time calculations marsden2. The fluid velocities are provided by the lattice Boltzmann solver described earlier employing fluid density and viscosity as g/cm and g/(scm), respectively, with a constant uniform velocity of cm/s imposed at the inlet and zero traction at the outlet (corresponding to ). A snapshot of the resulting steady-state solution in the domain is given in Figure 8 (right), where flow appears to strongly recirculate in the cavity (while remaining stagnant in the bottom cavity corners).

Using a spatial discretization of points (corresponding to cm) and a timestep of ms for the convection-diffusion PDE, the FC-based particle residence time calculations for cases with zero physical diffusion and non-zero physical diffusion ( cm/s) are presented in Figure 9. The simulation results demonstrate that the fluid gets trapped in the cavity and the stagnation corners as expected. In order to further validate the results of FC-based simulations with those provided by others marsden2, Table 1 additionally presents various measures/metrics computed over the entire domain as well as just the cavity. Defining as the effective length of time of a simulated “cycle” (of which there are ) and as the volume (area) of the region of interest (the full domain or only the cavity), two additional residence time measures and have been proposed marsden2 as

 RT1=1T∫NcyclesT(Ncycles−1)T1VΩi∫Ωiτ(x,t)dΩdtandRT2=VΩi1T∫NcyclesT(Ncycles−1)T∫Γiv(x,t)⋅ndΓdt, (30)

where the denominator in represents the average flow into the region of interest (e.g., across the boundary of the inlet or cavity opening where is the outward unit normal). Table 1 presents the corresponding values of these measures produced by FC-based simulations, which are well in agreement with those provided by others marsden2. The significance of these metrics has been discussed elsewhere marsden2.

Figure 10 presents snapshots at an s of an FC-based simulation of a dye injected across the inlet for s (or 2 cycles) for both zero physical diffusion and non-zero physical diffusion ( cm/s), where the latter illustrates early signs of the physical diffusion. Figure 11 presents corresponding snapshots at s, where the physically-diffusive configuration leads to larger concentrations of dye particles being trapped in the cavity (particularly the bottom-left stagnation corner, where a non-diffusive dye may take much longer to reach the area—if it does all).

### 4.2 Pulsatile flow through an aortic dissection (3D axisymmetric)

Calculation of particle residence times in a false lumen of an aortic dissection is physiologically significant since promoting thrombus (blood clot) formation in a false lumen (linked to PRTs) is of clinical interest erbel1993effect. The 3D axisymmetric computational domain of a simplified partially-grafted aortic dissection is illustrated in Figure 12, where both rigid walls and compliant walls (governed by the elastic PDEs of Equation (24)) are indicated. The fluid-structure immersed boundary-lattice Boltzmann solver described earlier is employed to simulate a physiological pulsatile flow and generate a flat velocity profile at the aortic inlet corresponding to a heart rate of bpm and a cardiac output of L/min. Fluid parameters correspond to and . Snapshots at various times of the pulsatile flow in the domain during a single cardiac cycle (of length s) is presented in Figure 13.

Using a spatial discretization of points (corresponding to cm) and a timestep of ms for the convection-diffusion PDE, the FC-based particle residence time calculation after many dozens of cardiac cycles for zero physical diffusion is presented in Figure 14. As expected, there are high residence times within the false lumen (between the septum and the aortic wall), particularly within the rigid grafted area.

Figure 15 presents snapshots of an FC-based simulation of a dye injected across a small part (the central 0.5 cm) of the inlet of the dissected region (within the true lumen) for s (or 5 cycles) for both zero physical diffusion and non-zero physical diffusion ( cm/s). The latter demonstrates how a diffusive dye (or drug) fills up much of the true lumen of the aorta. Figure 16 presents corresponding snapshots several cardiac cycles after the dye injection ceases, where in both cases the dye passes through the free end of the septum, remaining trapped even after several additional cycles (a similar effect is seen in the PRT calculations of Figure 14).

### 4.3 Non-Newtonian flow through a Fontan circulation (3D)

Fontan patients don’t have the right side of the heart, and blood flow is pushed from superior vena cava (SVC) and the inferior vena cava (IVC) to the lungs in a passive way via a cross-shaped graft (known as a “Fontan graft”) erbel1993effect. The 3D computational domain of a Fontan circulation is given in Figure 17 (left), where the top inflow is from the SVC and the bottom inflow is from the IVC. The two outlets flow towards the left and right lungs. The non-Newtonian fluid velocities are provided by the lattice Boltzmann solver described earlier with rigid walls and zero traction at the outlet (the shear rate viscosity curve of the shear-thinning Carreau–Yasuda model matching with Fontan patients is similar to wei2020). Blood fluid parameters correspond to , with a steady source corresponding to a cardiac output of 3 L/min. Further details of this configuration and its physiological and clinical relevancy are provided elsewhere wei2020. A snapshot (depth mid-plane cross-section) of the steady-state non-Newtonian flow velocity in the domain is presented in Figure 17 (right).

Using a spatial discretization of points (corresponding to cm) and a timestep of ms for the convection-diffusion PDE, the FC-based particle residence time calculation for zero physical diffusion is presented in Figure 18. The residence time is high at the area between the two outlets. Figure 19 presents snapshots during the initial cycles of an FC-based simulation of a dye injected across a subset of the top inlet for s (or 20 cycles) for both zero physical diffusion and non-zero physical diffusion ( cm/s). As expected, the physically-diffusive case begins to fill up the entire inlet of the Fontan graft. Figure 20 presents corresponding snapshots after several cycles, where in the physically-diffusive case, the dye particles diffuse into the area between the two outlets.

## 5 Conclusions

This work represents the successful construction of a new pseudo-spectral (FC-based) methodology for conducting dye simulations and non-discrete PRT calculations. Studies analyzing errors and convergence (with analytical solutions of both forced and unforced systems) attest to high-order accuracy without the incurrence of numerical diffusion errors (which are well-known to exist for commonly-used (low-order) finite element or finite difference-based methods kalinowska2007; sankar1998; mullen). Successful applications of the new solver to physically-realistic and physiologically-relevant case studies demonstrate the efficacy and versatility of the proposed approach. Although velocity data is provided in this work by numerical solutions (facilitated here by an FSI immersed boundary-lattice Boltzmann solver), the post-processing nature of the dye/PRT models can also apply to velocity data collected through suitable experimental techniques (e.g., particle image velocimetry) or clinical imaging modalities (e.g., cardiac 4D flow MRI or ultrasound image velocimetry). Analysis on the accuracy of the current method applied to such experimental data (e.g., compared with experimental dye tracing) is a subject of future work.