## Abstract

A finite-volume method has been developed that can deal accurately with complicated, curved boundaries for both two-dimensional and three-dimensional axisymmetric advection–diffusion systems. The motivation behind this is threefold. Firstly, the ability to model the correct geometry of a situation yields more accurate results. Secondly, smooth geometries eliminate corner singularities in the calculation of, for example, mechanical variables and thirdly, different geometries can be tested for experimental applications. An example illustrating each of these is given: fluid carrying a dye and rotating in an annulus, bone fracture healing in mice, and using vessels of different geometry in an ultracentrifuge.

## 1. Introduction

There are many physical and biological situations that can be modelled by the advection–diffusion equation:
1.1
where *q*(*x*,*y*,*z*,*t*) is the time (*t*)- and space (*x*,*y*,*z*)-dependent concentration of some chemical or density of some dispersed population, *κ*(*x*,*y*,*z*)≥0 is the capacity function, *D*(*x*,*y*,*z*,*t*) is the diffusion coefficient and **u**(*x*,*y*,*z*,*t*) is the velocity vector. Equation (1.1) is posed on a spatial (physical) domain or and for a time interval *I*=[*t*_{0},*t*_{end}]. The problem is complemented with initial and boundary conditions. In the case when , we assume that the problem is axisymmetric such that it can be reduced to an equivalent problem on its cross-section in .

Note that equation (1.1) is a simple example and the models normally encountered in nature, and later in this paper, are more complicated with multiple variables *q*=(*q*_{1},…,*q*_{l}) involved, upon which *D* and **u** can depend and with additional reaction terms between the variables. Equation (1.1) and its generalizations are so-called conservation laws, mimicking the conservative nature of the systems being modelled. For equation (1.1) in particular, integrals of over a spatial region are conserved in time up to fluxes through the boundary of that region. The finite-volume method is designed to maintain the conservation properties of such equations and hence is a natural approach for the computation of approximate solutions (e.g. LeVeque 2002).

Finite-volume methods can be employed on general spatial domains. Unstructured grids are good for boundary fitting of complex boundaries but are typically difficult to implement and can be computationally intensive. The need for fast numerical solution techniques together with a reasonably straightforward implementation of the computational scheme naturally leads to the use of a Cartesian grid on a rectangular spatial domain for equation (1.1). This latter approach is frequently used, although the physical situation might require a non-rectangular spatial domain *Ω*. Some groups improve on the allowed shape of *Ω* by using a non-intersecting union of these large rectangles to make up their domain (e.g. Geris *et al.* 2006). There are several disadvantages to this approach. Firstly, the rectangular domain may not faithfully represent the geometry of the situation modelled, thus creating reliability issues. The magnitude of this problem depends on the specific case and there is no real way of evaluating how large this effect might be. Secondly, if mechanical variables, such as stresses, are introduced into the mathematical model, then corner singularities can arise owing to the hard corners in the non-convex geometry that may arise from the union of rectangles. Mechanical forces play an important role in many biological processes including the example that will be considered here, namely bone fracture healing. The influence of the mechanics on the bone regeneration process will not be captured accurately by the model if singularities are present in the calculation of the mechanical stimuli. Thirdly, if a system is being designed for experimental or any other use, it would be desirable to investigate different geometries explicitly so that one which produced optimal results could be used.

The basis of the approach followed in this study is that of Cartesian cut cells in two dimensions developed by Calhoun & LeVeque (2000). The computational domain is a rectangle made up of grid cells (finite volumes) of dimensions *Δx*×*Δy*. The physical domain *Ω*, i.e. the domain on which equation (1.1) is to be solved, is a subset of the computational domain. There are uncut cells, i.e. cells that are either completely inside or completely outside *Ω*, and *cut cells*, i.e. cells that have been cut by a single straight line. By putting these cells together, an approximation to the physical domain *Ω* can be achieved figure 1. We will make use of the capacity function *κ* of the problem, extended to the whole computational domain, in order to encode the fraction of a grid cell that is part of *Ω*. A characteristic of this approach is that the numerical solution approximates the average of the solution in the Cartesian grid cells intersected with *Ω*. At times these averages are interpreted as solution approximations in the centroids of those intersections. The approach of Calhoun & LeVeque (2000) can be traced back to earlier work, e.g. Berger & LeVeque (1989), where the Euler equations in irregular two-dimensional spatial domains are considered. Even earlier, Cartesian cut-cell methods were applied in aerospace research for the solution of the transonic full potential equation in general domains (e.g. Wedan & South 1983). More current work in this direction is, for example, the *h*-box method for hyperbolic conservation laws (e.g. Berger *et al.* 2003).

There is a substantial body of literature dealing with cut-cell approaches to various applied problems. One slightly different line of research is that of Colella *et al*. (2006). Johansen & Colella (1998) use a cut-cell or ‘Cartesian-grid embedded boundary’ method to solve the (time-independent) Poisson equation with variable coefficients and Dirichlet boundary conditions on irregular-shaped domains *Ω* in one dimension and two dimensions. Their approach employs a cell-centred finite-volume discretization of the problem where the numerical solution approximates the exact solution in the centres of the Cartesian grid cells (and not in the centroids of those cells intersected with *Ω*). For certain geometries, this leads to cut cells where the solution is approximated outside *Ω*. This could be seen as a disadvantage but on the other hand it is key to arriving at linear equation systems whose condition is relatively unaffected by small cut cells. The Dirichlet boundary data are incorporated using a quadratic polynomial to approximate the solution and lead to one-sided finite differences for the computation of the boundary fluxes. The scheme displays second-order accuracy for test examples. The approach above has been extended to the time-dependent case (McCorquodale *et al.* 2001), and the heat equation for irregular two-dimensional spatial domains is considered. Here, a dedicated implicit, second-order time-stepping scheme (without time-step size control) is employed. For a particular parameter choice, this scheme encompasses the Crank–Nicolson method. The method is extended to three dimensions in Schwartz *et al.* (2006). A similar approach is taken by Colella *et al.* (2006) for hyperbolic conservation laws. The spatial discretization uses a flux-redistribution approach to ensure an overall conservative scheme despite the use of some non-conservative algorithmic ingredients.

In the case where the irregular domain *Ω* is rather simple, e.g. a circle or a ring or similar, then logically Cartesian grids can be employed, often more easily (for a current review see Calhoun *et al.* 2008).

This study modifies the two-dimensional Cartesian cut-cell technique of Calhoun & LeVeque (2000) for use in a method of lines (MOL) approach. It is extended to allow for Dirichlet boundary conditions. The use of limiter functions yields accurate advective flux approximations and avoids the introduction of oscillations into the numerical solution. Finally, the spatially two-dimensional approach is extended to be applicable for three-dimensional axisymmetric systems. Whereas Calhoun & LeVeque (2000) consider the discretization of diffusion and advection individually and combine them in a fractional-step time-stepping scheme, the MOL approach (e.g. Hundsdorfer & Verwer 2003; Gerisch & Chaplain 2006) treats the discretization of both parts in the same framework resulting in a single system of ordinary differential equations (ODEs), the MOL–ODE. This has the advantage that the MOL–ODE can be, and is, solved with advanced time-stepping schemes of high order and with automatic time-step control to ensure an efficient and accurate solution.

## 2. Derivation of an MOL cut-cell discretization

In this section the derivation of the cut-cell discretization is presented for the model problem equation (1.1) for two spatial dimensions (*x*,*y*). As already described in the introduction, the method of Cartesian cut cells involves defining a grid that covers a computational domain greater than the physical domain *Ω*. The grid cells are denoted by *Ω*_{ij}; *i* runs in the *x*-direction and *j* in the *y*-direction. With (*x*_{i},*y*_{j}) we denote the centre of grid cell *Ω*_{ij}. When the grid is placed over *Ω*, there will be cells *Ω*_{ij} that lie totally outside, those that are totally inside and those that are cut by part of the boundary of the physical domain. Following Calhoun & LeVeque (2000), we assume that cut cells are cut with a single straight line and extend the capacity function *κ* of equation (1.1) from *Ω* to the whole computational domain by defining it to be zero outside *Ω*. Then for each cell, we can define its average capacity as
It readily follows that *κ*_{ij}=0 for all cells *Ω*_{ij} totally outside *Ω*. In many applications, we have *κ*≡1 in *Ω* and then *κ*_{ij}=1 for *Ω*_{ij} totally inside *Ω* and 0<*κ*_{ij}<1 for cut cells. Hence, in that case *κ*_{ij} is the fraction of *Ω*_{ij}, which is part of the physical domain. The quantities that we wish to compute are time-dependent approximations to the spatial average of the solution *q*, weighted with the capacity function, in each grid cell, i.e.
If *κ*≡1 in *Ω*, then *Q*_{ij}(*t*) simply approximates the average of *q* at time *t* in the intersection . The set-up of a system of ODEs for the averages *Q*_{ij}(*t*), based on the partial differential equation (1.1) at hand, constitutes the first step of the MOL. The initial values *Q*_{ij}(0) for this MOL–ODE are taken from the initial data that accompany equation (1.1). In this paper, we are primarily concerned with the set-up of the MOL–ODE. For the second step of the MOL, the time integration of the MOL–ODE, we employ the linearly implicit Runge–Kutta method of order four ROWMAP (Weiner *et al.* 1997). This fourth-order method employs Krylov techniques for the solution of the large linear equation systems arising in its stages and is well suited for the numerical solution of large and stiff ODE systems such as our MOL–ODE.

For a grid cell *Ω*_{ij} with zero capacity *κ*_{ij}, as is the case in cells completely outside the physical domain, we define the corresponding ODE and initial condition
2.1a
For all other grid cells, the ODE is obtained by accounting for the fluxes through all the grid-cell boundaries. Following the finite-volume methodology in Calhoun & LeVeque (2000), this then amounts to
2.1b
In this equation, *F*_{i+1/2,j},−*F*_{i−1/2,j},*F*_{i,j+1/2} and −*F*_{i,j−1/2} are the average (per unit length) fluxes across the right, left, upper and lower edges of *Ω*_{ij} in the outward normal direction. This means, for example, that *ΔyF*_{i+1/2,j} is the total flux out of *Ω*_{ij} through the right boundary edge. Conservation of mass dictates that *F*_{i−1/2,j}=*F*_{(i−1)+1/2,j} and *F*_{i,j−1/2}=*F*_{i,(j−1)+1/2} and consequently for each grid cell only the fluxes out across the right and upper edges must be computed. In cut cells, there is another boundary segment, namely that which coincides with or represents a part of the physical boundary ∂*Ω*. We denote the length of this segment by *l*_{ij} and the average (per unit length) flux out through this segment by *H*_{ij}. It should be noted that in the case of no-flux boundary conditions, *H*_{ij}=0 holds. Reaction terms can also be added to this ODE (e.g. Gerisch & Chaplain 2006).

In the following §2*a*,*b*, we treat the definition of the average fluxes through the edges of grid cells. These fluxes are in the direction for a right grid-cell edge, for an upper grid-cell edge and in the direction of the unit outward normal vector of *Ω* for the additional edge of length *l*_{ij} of a cut cell. The average flux *F* through an edge in the direction is the sum of an average diffusive, *F*^{D}, and an average advective, *F*^{A}, flux in the direction . The diffusive and the advective flux in direction in a point (*x*,*y*,*t*) in equation (1.1) are
respectively. The cut-cell approach requires the distinction between uncut- and cut-cell edges of grid cells. For the diffusive flux computation, an edge is an uncut-cell edge if and only if the two adjacent grid cells are uncut grid cells, otherwise that edge is a cut-cell edge for diffusion. For the advective flux computation, an edge is an uncut-cell edge if and only if the two adjacent grid cells upstream and the one downstream are uncut grid cells, otherwise that edge is a cut edge for advection (see figure S1 in the electronic supplementary material for an illustrative sketch). Subsequently, the implementation of boundary conditions for equation (1.1) (§2*c*), and an extension of the two-dimensional cut-cell approach for the three-dimensional axisymmetric problems (§2*d*), are presented.

### (a) Flux computation on uncut-cell edges

In the case of an uncut-cell edge, on either side of that edge there is an uncut (full) grid cell *Ω*_{ij} (in the case of the advective flux computation, there are even two uncut grid cells upstream) with the average concentration of *q* given by *Q*_{ij}. We can also interpret *Q*_{ij} as an approximation of *q* in the grid-cell centre (*x*_{i},*y*_{j}), which is equal to the centre of the mass of *Ω*_{ij}. For these reasons, we therefore apply standard finite-volume approximation techniques to define the fluxes for the MOL–ODE.

#### (i) The diffusive flux on uncut-cell edges

For a right edge, the outward normal vector is and the diffusive flux through that edge is −*Dq*_{x}. We follow Calhoun & LeVeque (2000) and approximate the derivative *q*_{x} by central differencing across the edge, then evaluate *D* in the edge centre and consequently define
2.2
Accordingly, for an upper edge, the normal is given by , and we define
2.3

#### (ii) The advective flux on uncut-cell edges

In many systems, steep gradients exist in the exact solution and within the numerical solution process, these cause oscillations to appear in the numerical solution if the advection term is not carefully discretized. These oscillations may result in negative numerical solution values, although the exact solution is non-negative, for example modelling a concentration. These negative values create major problems in many systems, particularly those which also involve reaction terms that may render the differential equation unstable if they are evaluated for negative (numerical) solution values and hence amplify the oscillatory behaviour. To combat this, higher-order upwinding discretizations together with so-called limiter functions have been introduced (see Hundsdorfer & Verwer 2003 for more details). In the following approach, we make use of the van Leer (1974) limiter. The introduction of limiters effectively means that a locally appropriate combination of higher- and first-order upwind difference schemes is used to create an accurate but non-oscillatory solution.

Let us firstly consider the advective flux through a right uncut-cell edge, . The definition of that flux depends on the sign of the normal velocity in the centre of the uncut-cell edge, **u**(*x*_{i}+*Δx*/2,*y*_{j},*t*)⋅(1,0)^{T}=*u*_{1}(*x*_{i}+*Δx*/2,*y*_{j},*t*). For a positive velocity, *u*_{1}(*x*_{i}+*Δx*/2,*y*_{j},*t*)>0, we use
2.4aand
2.4b
together with the van Leer limiter function
2.5
Here *r*_{i+1/2,j} is a measure of the monotonicity of the grid data. In particular, if *Q*_{ij} is an extrema with respect to its neighbours *Q*_{i−1,j} and *Q*_{i+1,j} then *r*_{i+1/2,j}=0 and also *L*(*r*_{i+1/2,j})=0. This implies that in this case the flux is computed using the first-order upwind discretization. If *L*(*r*_{i+1/2,j})≠0 then a higher-order discretization is employed on the edge. From the above equations, it is also apparent that the two grid-cell values *Q*_{i−1,j} and *Q*_{ij} upstream and the grid-cell value *Q*_{i+1,j} downstream of the edge are used in the definition of . These grid cells must therefore be uncut grid cells completely within the physical domain *Ω*. If they are not, then the edge under consideration is a cut-cell edge (for advection) and a flux computation appropriate for this case (see further below) must be employed.

For a negative velocity, *u*_{1}(*x*_{i}+*Δx*/2,*y*_{j},*t*)<0, the upwind direction changes but the equations for the advective flux are similar and given by
2.6aand
2.6b
Here again, two values, *Q*_{i+1,j} and *Q*_{i+2,j} upstream and one value, *Q*_{ij}, downstream of the edge are used in the definition of . It is straightforward to write down the definition of the advective flux on an upper uncut-cell edge, although we abstain from it here. A thorough investigation of the advective flux discretization as given above can be found in Gerisch (2001).

There is a slight error in this approach caused by the splitting of the flux in the *x*- and *y*-direction. This is illustrated in figure 2, in which we assume a constant advection velocity **u**=[*u*_{1},*u*_{2}] with *u*_{1},*u*_{2}>0 and assume that the concentration *q* is constant in each grid cell *Ω*_{ij} at the value *Q*_{ij}. Then, over a small time span *δt*, the differences between the actual and computed (first-order upwind) fluxes in and out of *Ω*_{ij} as well as their sum are given by
These errors decrease with *δt* and so smaller time steps (resulting in longer integration times) will produce more accurate results such that the temporal integration can deal with these errors to the required degree of accuracy. If the concentration gradient is small, then all the concentration terms will almost entirely cancel each other out, thus relaxing the constraints on *δt* and reducing the integration time.

### (b) Flux computation on cut-cell edges

The situation is now different from that for uncut-cell edges: there llwi be at least one adjacent (or next to adjacent for the advective flux) grid cell, which is a cut cell. For such a cut cell *Ω*_{ij}, the value of *Q*_{ij} approximates the average of *q* in . So we should not interpret *Q*_{ij} as an approximation of *q* in the grid-cell centre (*x*_{i},*y*_{j}) but as an approximation in the centre of mass of .

If we consider the average flux *F*_{i+1/2,j} on a right cut-cell edge, then *F*_{i+1/2,j} must be defined such that *ΔyF*_{i+1/2,j} is the total flux through that edge. This total flux accumulates only on that part of the grid-cell edge which is within the physical domain. If we denote the length of this part by *h*_{i+1/2,j}, then as in Calhoun & LeVeque (2000), a consistent definition is obtained by letting
2.7
where is the average (per unit length) flux out of *Ω*_{ij} through the portion of length *h*_{i+1/2,j} of its right edge. Similarly, we define the average flux *F*_{i,j+1/2} on an upper cut-cell edge by
2.8

#### (i) The diffusive flux on cut-cell edges

For the diffusive flux, the same approach as that by Calhoun & LeVeque (2000) is taken where a bilinear interpolant
2.9
is fitted to the data *Q*_{ij} in the centres of the mass of four neighbouring cells of the cut-cell edge under question. Then can be differentiated analytically and evaluated at the centre (*x*_{e},*y*_{e}) of the part of the cut-cell edge within *Ω*. The diffusion coefficient is also evaluated at this point. In summary, we then obtain for a right cut-cell edge
2.10
and for an upper cut-cell edge with a different point (*x*_{e},*y*_{e}),
2.11

We note that if the approach outlined here for cut-cell edges is applied to uncut-cell edges, then the same flux value on the uncut-cell edge as before is obtained.

#### (ii) The advective flux on cut-cell edges

Here, we firstly consider the case that the second grid cell upstream of the cut-cell edge under consideration is outside *Ω*; this is a typical scenario near an inflow boundary part of *Ω*. Secondly, we consider the case that this grid cell is, at least partially, contained in *Ω*, a scenario often encountered near an outflow boundary part of *Ω*. Finally, the special procedure to compute the total cell flux owing to advection for cut cells *Ω*_{ij} with very small capacity 0<*κ*_{ij}≪1 is described. We add this last special case in order to obtain a discretization algorithm that can deal reliably with very small cut cells.

If the second grid cell upstream of the cut-cell edge under consideration is outside *Ω*, then it is not possible to apply the discretization with van Leer limiter owing to unavailable values of *Q*_{ij}. Hence, we resort to the first-order upwind discretization in that case and define for a right cut-cell edge
Here again (*x*_{e},*y*_{e}) denotes the centre of the part of the cut-cell edge within *Ω*. The value on an upper cut-cell edge is defined accordingly.

Also, if the second grid cell upstream of the cut-cell edge under consideration is, at least partially, inside *Ω*, then we have the available necessary values *Q*_{ij} to apply the discretization with van Leer limiter which we have used to define the advective flux on uncut-cell edges. However, those *Q*_{ij} should be considered approximations to *q* in the centres of mass of their respective *Ω*_{ij}. Nevertheless, we apply the equations for uncut-cell edges but with the velocity *u*_{1} on a right cut-cell edge (respectively, *u*_{2} on an upper cut-cell edge) evaluated in (*x*_{e},*y*_{e},*t*), where (*x*_{e},*y*_{e}) again denotes the centre of the part of the cut-cell edge within *Ω*. In particular, we define for a right cut-cell edge and the case *u*_{1}(*x*_{e},*y*_{e},*t*)>0, cf. equation (2.4),
2.12
with *r*_{i+1/2,j} and *L*(*r*) defined as before. for *u*_{1}(*x*_{e},*y*_{e},*t*)<0 is defined following equation (2.6) and the value on an upper cut-cell edge is defined accordingly. This treatment of the advective flux on these particular edges will in general reduce the order of accuracy to one locally (similar to when the limiter function reduces the accuracy to first order in the vicinity of local extrema of the solution).

The advective fluxes computed for a grid cell *Ω*_{ij} are added up and divided by *κ*_{ij} (see equation (2.1b)). Whereas this division does not cause a difficulty in the case of uncut grid cells, it may pose a source of error for small cut cells where 0<*κ*_{ij}≪1. Therefore, for these small cut cells, we require a different approach since otherwise the discretization may either lead to numerical instability or make it necessary to limit the time step in the time integration to a very small value.

In the following, we consider for each *Ω*_{ij} that part of the right-hand side of equation (2.1b) which arises from the discretization of advection and denote it by ,
This value represents the average (per unit area) change of mass of *Q* in grid cell *Ω*_{ij} owing to advection. Our basic assumption for cut cells *Ω*_{ij} with 0<*κ*_{ij}≪1 is that will be similar to the value in adjacent grid cells *Ω*_{lm} as the distance between their respective centres of mass, and , will be small. Therefore, we fit a surface through the (known) values in the adjacent cells and evaluate this surface at the centre of mass of *Ω*_{ij} to find an approximate value of . This may cause some minor problems with the conservation of mass but owing to the fact that this approach is only taken for cells with low values of *κ*_{ij} and the mass within a cell is *κ*_{ij}*ΔxΔyQ*_{ij}, this slight variation will be minimal (see the electronic supplementary material, Section 2). Also this modification has been shown to have a negligible effect on the numerical errors in test cases but it results in an often substantial reduction of CPU time for the time integration of the resulting MOL-ODE (see also Section 2 in the electronic supplementary material). If in later applications the lack of conservation of the scheme is deemed too large then a flux-redistribution approach such as that in Colella *et al.* (2006) could be implemented. Numerical tests (not presented in detail in this paper) showed that a suitable condition for a cut cell *Ω*_{ij} to be counted as small was that *κ*_{ij}<0.01.

In order to obtain an accurate approximation of the unknown value in a small cut cell, the surface should be fitted to the (known) values of in as many as possible adjacent (cut and non-cut) cells. To this end, a procedure with at most five steps is applied to find for all small cut cells by interpolation: in step *k*, we determine all still unknown values where the are already known in at least #_{k} adjacent cells of the small cut cell *Ω*_{ij}. The pairs (*k*,#_{k}) are given by (1,4), (2,3), (3,3), (4,2), (5,1) so that in the end all have been determined.

### (c) Implementation of boundary conditions

Boundary conditions for problem (1.1) can apply to cut or uncut upper, lower, left or right edges of a grid cell *Ω*_{ij} or to the special single straight edge that may make *Ω*_{ij} a cut cell. The effect of these *boundary edges* is captured in equation (2.1b) by the term (one for each boundary edge of *Ω*_{ij})
The numerical treatment of all those boundary edges follows the same strategy. We consider one such edge of length *l*_{ij} and centre (*x*_{e},*y*_{e}) and explain the definition of the average flux *H*_{ij} out of *Ω*_{ij} through that edge.

If the flux on the boundary edge is prescribed as part of problem (1.1), then we define *H*_{ij} simply as the prescribed value in (*x*_{e},*y*_{e}). In particular, no-flux boundary conditions lead to *H*_{ij}=0.

For Dirichlet boundary conditions, the value of *q* is prescribed on the boundary edge and we make use of the value in its centre, *Q*_{D}:=*q*(*x*_{e},*y*_{e},*t*). Furthermore, we let be the unit outward normal vector of the boundary edge. The average advective flux through the boundary edge is defined by first-order upwinding
The situation is slightly more complicated for the average diffusive flux. In that case, we assume that this is aligned with the boundary edge but outside *Ω* with a square ‘ghost’ grid cell with side of length in which *q* takes the constant Dirichlet value *Q*_{D}. The centre of the ghost cell is found by first projecting the centre of mass of *Ω*_{ij} onto the boundary edge extended to a straight line, defining a distance vector and then continuing in this direction for a distance *δ*/2. Hence, the distance between the centre of mass of the ghost cell and that of *Ω*_{ij} is given by *α*+*δ*/2. The average diffusive flux through the boundary edge is then defined by central differencing across the boundary edge

### (d) Extension to three-dimensional axisymmetric problems

Thus far the numerical technique has been developed for problem (1.1) in two spatial dimensions. A considerable number of three-dimensional problems, including two of those treated in §4, have the property of axisymmetry and can equivalently be reduced to a two-dimensional problem on corresponding to half of their cross-section. For instance, a three-dimensional axisymmetric problem (1.1) reduces to
Here *x* is the axial coordinate and *r*≡*y* the radial coordinate. It should be noted that this system is still in conservation form but in particular *q* is multiplied by the radial coordinate *r* on the left-hand side. By introducing ‘radius-weighted’ averages *Q*_{ij} of the solution and following the approach in Gerisch & Geris (2007), we can apply the numerical techniques developed above to this modified problem.

## 3. Numerical tests of the cut-cell discretization

In this section, particular aspects related to the cut-cell discretization described in §2 are discussed. Convergence studies were carried out to test the influence of the various constituents of the model problem (1.1) as well as the relative size of the cut cells on the spatial discretization. Furthermore, we investigated numerically the influence of the introduction of the van Leer limiter in the spatial discretization of the advective term of equation (1.1) on the overall accuracy and convergence speed. The time integration of the MOL–ODE system was carried out with ROWMAP (Weiner *et al.* 1997) using a sufficiently strict tolerance of 10^{−8} such that errors of the time integration process could be considered negligible compared with the spatial discretization errors that are of interest for the convergence studies in this section.

### (a) Convergence of the cut-cell discretization

Convergence tests were carried out for the numerical solution of problem (1.1) with (i) diffusion alone (ii) advection alone and (iii) diffusion and advection together. A simple model of a pulse of suspension moving along a two-dimensional pipe of width 1 m and at an angle *α*=*π*/4 to the positive *x*-axis with constant velocity *v* and/or dispersing with constant diffusion coefficient *D* was used as the basis. Consequently, the velocity vector of equation (1.1) is given by . The initial shape of the pulse is captured in the initial condition
(figure 3*a*). Here, *d*(*x*,*y*) is the distance of (*x*,*y*) from the line perpendicular to the pipe and through the centre of the computational domain. We assume zero flux boundary conditions along the pipe. The physical domain *Ω* is embedded in a computational domain of size of 5 m by 5 m (figure 3*a*) and simulations were run for 8 s. The spatial discretization is also supposed to work well for small cut cells in the grid. To test the behaviour of the discretization in such cases, the lower edge of the pipe was moved to three different positions to give cut cells of half capacity, one eighth capacity and one eight hundredth capacity as shown in figure 3(*b*).

Figures 4 and 5 show the results of the convergence study. We use grids with square cells of width and consider the error in space at the final time. To calculate this error for the different cell widths *Δx*, a reference solution was needed. For the pure advection case, the analytical solution of the model problem, i.e. the appropriately displaced initial condition of the problem (figure 3*a*) was used as the reference solution. For the case involving diffusion alone, the results of a computation on the grid with *Δx*=5/150 were taken to be sufficiently accurate and used as the reference solution. Finally, the reference solution for the combined case of diffusion and advection was obtained by analytically displacing the reference solution of the pure diffusion case at the final time to account for the advection. To calculate the reference solution for a grid with cell width *Δx*>5/150, the concentrations at the centre of mass of each cell were linearly interpolated from the finer reference grid. As a measure of error, we approximate the volume weighted *L*^{2}(*Ω*)-norm of the error, function *e*(*x*) from the difference between numerical, *Q*_{ij}, and reference solution, , values on the grid:
3.1

From figure 4*a* it is obvious that the algorithm applied to problem (1.1) with diffusion only converges almost with order two for decreasing spatial grid width. Furthermore, no significant dependence of the error on the size of cut cells is observed. This graph gives the least errors in all the tests shown in figures 4 and 5 and so we checked that the level of temporal tolerance (10^{−8}) for ROWMAP is sufficiently strict by recalculating the solution on the 140×140 grid with a stricter temporal tolerance of 10^{−10}. Both errors were found to deviate by less than 3×10^{−5}, thus showing that the original temporal tolerance (10^{−8}) was strict enough.

The algorithm applied to problem (1.1) with advection only (figure 4*b*) converges at a slower rate and when it is applied to the combined diffusion–advection problem, figure 5*a*, the error lies between the error in solely the diffusion case and that in solely the advection case. This is due to the fact that a present diffusion term smoothes the solution (i.e. decreases the gradients). This means that the discretization of the advection (if present) can make more use of the higher-order unlimited approximation of the advective flux than in cases with less or no diffusion where steep gradients together with the limiter function drive the advective flux approximation closer to a first-order upwind discretization; this effectively adds *numerical diffusion* to the scheme in order to avoid having oscillations occur in the numerical solution. There is also a slight tendency for very small cells to have larger errors, which is the reason for the slight increases in the error at various grid sizes since the upper boundary cuts the cells in different ways. From figure 5*b*, it can be seen that the worst errors come from systems with large amounts of advection and little diffusion.

### (b) The effect of the van Leer limiter

One of the major differences between the method used in Calhoun & LeVeque (2000) and our method is that we treat diffusion and advection in the same framework (MOL) resulting in a single ODE system after spatial discretization. This allows the use of off-the-shelf, high-order time integration schemes such as ROWMAP. Calhoun & LeVeque (2000) treat diffusion and advection independently and combine both using a fractional-step method. This usually limits the temporal accuracy to order two. For advection, they employ the cut-cell approach with the Clawpack software that provides high-resolution discretization of hyperbolic conservation laws such as, for example, problem (1.1) without diffusion. To also have a better than first-order upwind discretization of the advection part in our combined approach to the advection–diffusion problem (1.1), we make use of a van Leer limited, higher-order upwind discretization of the advective flux. This is an important improvement over the simple upwinding but it makes the calculation process slower by up to a factor of two on the same grid. Therefore, the advection convergence test with half-cut cells was also run for two different velocities with the van Leer limiter function replaced by *L*(*r*)=0, leading to the first-order upwind discretization of the advection term, to show its effect. The results are shown in figure 5*b*. It can be seen that there is a great improvement in accuracy when using the discretization with van Leer limiter: convergence faster than first order and much lower errors are observed. It is therefore well worth spending the additional computational effort to compute the van Leer limited advection discretization as otherwise with first-order upwinding very fine grids would be required to obtain equally accurate results but these would take much longer to compute.

## 4. Three example cases

This section contains examples that demonstrate the added value of different aspects of the approach presented in this paper. The example of a fluid carrying a dye demonstrates the influence of geometrical approximations on the computed results. The second example demonstrates how the use of unsmooth geometries might induce corner singularities that could adversely influence the model’s predictions. Finally, the added value of the approach in this study on the design of optimal vessel geometries for ultracentrifuge experiments is shown.

### (a) Fluid rotating in an annulus

#### (i) Problem

This example illustrates how approximating the geometry of a situation can greatly affect the results. We firstly consider a shallow circular annulus that is rotating slowly so that a viscous liquid contained within it will rotate with constant angular velocity. If a dye is then added to the liquid at a particular time, it will diffuse through the liquid while also being carried round with an advection term owing to the motion of the liquid that it is suspended in. Obviously, this situation could be simplified by changing frames but we remain in our current frame for the purpose of demonstration, as in more complicated cases that approach may not be possible. Here, we consider two possible geometries : firstly, *Ω* is obtained by approximating the annulus by a square structure and secondly, as favoured in the cut-cell approach of this paper, using cut cells to approximate the circular shape of the annulus. In both cases, we use no-flux boundary conditions on the boundary of *Ω*.

#### (ii) Results

Figure 6 contains the results for both the square (top row) and the cut-cell (bottom row) geometry for an annulus of outer radius *L* and inner radius 0.25*L*, rotating at 0.1 r.p.m. anticlockwise and containing a dye with diffusion constant 0.01*L*^{2} min^{−1}. In both cases, square grid cells with *Δx*=*Δy*=0.02*L* were used. The initial conditions are given by
where *θ*(*x*,*y*) is the angle of (*x*,*y*) measured round the annulus. They are shown in the left column of figure 6 and represent the situation that the dye (dark colour) is slightly diffused in the viscous fluid (light colour) already. It can clearly be seen that these give very different results. In the square approximation, the dye is pushed out radially as it rounds the edges of the inner square creating the more concentrated area near the outer edge. Also, dye diffuses out into the corners of the larger square and hence the total amount of dye within the actual annulus appears to decrease over time. This is obviously not physical as it violates mass conservation.

### (b) Bone fracture healing in mice

#### (i) Problem

The process of bone fracture healing is a complicated one involving many types of cells and biological material. Several different models have been created for this process, normally in the form of a set of differential equations with chemical reaction terms, diffusion terms and advection (taxis, active transport) terms. A typical fracture is characterized by the formation of a callus. This is a loose tissue type that is formed from the blood clot that emerges after fracture and contains most of the cells and factors that are needed for a successful fracture repair. It is within this callus that all the reactions, diffusion and active transport occur. In previous computer simulations of bone fracture healing, this volume was constructed crudely from rectangular sections as shown in figure 7*c* (see Geris *et al.* 2006). The method described in this paper is able to more accurately model the geometry of the situation (figure 7*b*).

For this example, the model equations established by Bailón-Plaza & van der Meulen (2001) are applied to simulate the healing process of a semi-stabilized murine tibia fracture in an axisymmetric callus geometry (cf. Gerisch & Geris 2007). We refer to Geris *et al.* (2006) for more details, initial and boundary conditions and parameter values. In this model, the form of the active transport (taxis) velocity of variable *q*_{j} is **u**_{j}=∇*g*, where *g* is a linear combination of the other variables *q*_{i}, *i*≠*j*. To approximate ∇*g* on the edge of a cell, exactly the same formulation as that used to determine ∇*q* for the diffusive flux is used. As indicated above, we consider here two different geometries for the callus, one of which is very similar to that of Geris *et al.* (2006) (figure 7*c*) and in the other one, the curved callus boundary is approximated more realistically using cut cells (figure 7*b*).

#### (ii) Results

Figure 8 shows the time evolution of the density of bone extracellular matrix for the two different callus geometries. (A similar figure but giving the stem cell density is provide in the electronic supplementary material.) It can be seen that the results are very similar thus justifying the approach of Geris *et al.* (2006). Slight differences are seen in the total healing time, which can be attributed mainly to the difference in the modelled callus volume.

The more realistic geometric model (figure 7*b*), of the fracture callus allows for an improved incorporation of mechanical loading in the mathematical model. Typically, mechanical stimuli such as tissue stresses are involved in the determination of cell fates and in the production of tissues and growth factors by these cells. Incorrect calculations of these mechanical stimuli will lead to incorrect predictions of the biological events by the mathematical model. Using geometries with re-entrant corners may result in singularities in the calculations of these mechanical stimuli and furthermore, the global mechanical behaviour of a structure depends on its actual shape. The presence of additional mass in certain places will enhance the overall stiffness of the structure. With the cut-cell approach, the number of re-entrant corners in the physical domain can be reduced (see the two domains in figure 7*b*,*c* considered here). Figure 9 shows the intermediate principal strain at time *t*=0 for both geometries under normal loading conditions. The calculated differences, both owing to corner singularities (arrow) and differences in the structure’s global stiffness, run up to 200 per cent, which would result in significant differences in the predicted regeneration process.

### (c) Different geometries in an ultracentrifuge

#### (i) Problem

The ultracentrifuge is a centrifuge optimized for spinning a rotor at very high speeds and capable of generating acceleration as high as 1 000 000*g* (9800 km s^{−2}). In an analytical ultracentrifuge, a sample being spun can be monitored in real time through an optical detection system to observe the evolution of the sample concentration versus the axis of rotation profile as a result of the applied centrifugal field. These observations are electronically digitized and stored for further mathematical analysis. The kinds of information that can be obtained from an analytical ultracentrifuge include the gross shape of macromolecules, the conformational changes in macromolecules and size distributions of macromolecular samples. Svedberg & Fåhraeus (1926) provided one of the initial estimates of the molecular weight of haemoglobin using the sedimentation rate and more recently this value has been used to determine the size of proteins in solution (see Selivanova *et al.* 2003).

In this example, we consider whether different geometries of the centrifuge tube could lead to more accurate results in studies of the type mentioned above. Spinning in a centrifuge creates a high, outward, artificial gravitational field. Any suspended particles are therefore accelerated outwards by this field but quickly reach their terminal velocity owing to the drag of the solution. The Mason–Weaver equation
describes the time evolution of the process of build-up of suspension towards the bottom of the tube. It is clear that this is advection–diffusion equation (1.1) with a diffusion constant *D* and a velocity field **u**=(0,0,−*sg*)^{T}, where *g* is the gravitational field and *s* is the sedimentation coefficient. In the simulations, we take *D*=1.94×10^{−9} m^{2} s^{−1} and *s*=2.6×10^{−13} s (roughly equivalent to a low concentration solution of haemoglobin—for high concentrations the variation of *D* with concentration becomes important). Using a centrifuge that produces a gravitational field of 81 000*g* enables the comparison of the results that would be expected from two different three-dimensional axisymmetric geometries of centrifuge tube. Geometry A is simply a cylinder of radius 0.5 cm and length 3 cm. Geometry B is a cylinder of length 2 cm and radius 0.5 cm with a cone of height 1 cm and a base of radius 0.5 cm on the end. The cross-sections of both are shown in the outermost plots of figure 10. Owing to the symmetries of these shapes only half of the cross-section needs to be considered. These tubes have large radii in comparison to their length so that the gravitational field *g* is constant across them. For the computational grid, we use grid widths *Δx*=*Δy*=1/20 cm.

#### (ii) Results

The middle plots of figure 10 show the central concentration variation at half-hourly intervals for the first 3 h for an initially uniformly distributed concentration, *q*(*x*,*y*,0)=1. There are slight differences between the two different geometries. Using an ultracentrifuge for a suspension with a known diffusion constant *D*, the experimental data can be fitted to the results of the simulation to find a value of the sedimentation coefficient *s*. As different geometries generate different concentration profiles, it is important to accurately capture the geometry of the centrifuge tube in order to make correct predictions of the sedimentation coefficient.

A simple method of determining the sedimentation coefficient *s* is to measure the absorption coefficient and then to fit this to the results of a numerical simulation. If the suspension is optically thin then the absorption coefficient is proportional to the integral across the tube diameter of the concentration *q*. To determine the sensitivity of this method, the simulation can be run with a slightly different value of the sedimentation coefficient, *s*_{2}:=1.01*s*_{1}, then observing the difference in the integrated concentration *q*. Figure 11 shows this difference for the geometries A and B. The integral of the absolute value of this difference for geometry B is greater than that for geometry A, hence the error in an experiment using geometry B will be smaller. By carrying out these simulations for a variety of geometries prior to the actual experiment, the centrifuge tube geometry can be designed so that the experimental results will be most accurate.

## 5. Conclusion

In this paper, a method of numerically solving the advection–diffusion equation in two-dimensional and three-dimensional axisymmetric space has been developed. An approach using rectangular grid cells that may be cut by a single straight line was used to deal with curved boundaries (cut-cell approach). There are three main advantages of this approach: more accurate results can be obtained, the possibility of singularities being obtained in the calculation of, for example, mechanical variables resulting from geometries with re-entrant corners is minimized and different experimental methods can be simulated to find those that will give the most accurate results. Examples of all of these have been given in the cases of fluid rotating in an annulus, bone fracture healing in mice, and the determination of the sedimentation coefficient in a centrifuge. These suitably demonstrated the usefulness of this approach. Computation of this algorithm was also reasonably quick with the most complicated example being the bone fracture healing that involved seven coupled variables being completed in only a few hours on a standard Pentium 4 processor.

We have presented a cut-cell approach within an MOL framework. One of our main motivations for this is that this framework allows us to use high-order time-stepping schemes (in our case ROWMAP, order 4), which employ temporal error estimation to adjust the time-step size automatically and which can consequently be expected to work efficiently towards the user-defined temporal accuracy goal. In comparison, fractional-step methods employed for diffusion–reaction problems, e.g. Calhoun & LeVeque (2000), are at most accurate to second order in time and hence with ROWMAP as a time integration scheme, we can expect a more efficient solution process for higher accuracy demands (and sufficiently smooth problems).

We note that the software created for the simulations in this work was conceived in a generic way. This allows the possibility of future extensions of the code to include locally refined grids and to go from two-dimensional to three-dimensional space.

## Acknowledgements

The authors would like to acknowledge the valuable and constructive comments from the referees. J.M.A.A., L.G., A.G. and C.J.S.Y. would like to thank the Nuffield Foundation for the award of an Undergraduate Research Bursary (Grant Ref: URB/35427) that supported this research. L.G. is a post-doctoral research fellow of the Research Foundation, Flanders (FWO Vlaanderen).

## Footnotes

- Received October 5, 2009.
- Accepted December 4, 2009.

- © 2010 The Royal Society