A numerical method for the multiphase viscous flow equations

Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A numerical method for the multiphase viscous flow equations James M. Osborne, Jonathan P. Whiteley ⇑ Computing Laboratory, University of Oxford, United Kingdom

a r t i c l e

i n f o

Article history: Received 22 December 2009 Received in revised form 3 June 2010 Accepted 20 July 2010 Available online 21 August 2010 Keywords: Multiphase viscous flow Finite element method Tissue development

a b s t r a c t A numerical technique is developed for the solution of the equations that govern multiphase viscous flow. We demonstrate that the equations can be written as a coupled system of Partial Differential Equations (PDEs) comprising: (i) first order hyperbolic PDEs for the volume fraction of each phase; (ii) a generalised Stokes flow for the velocity of each phase; and (iii) elliptic PDEs for the concentration of nutrients and messengers. Furthermore, the computational domain may vary with time for some applications. Appropriate numerical methods are identified for each of these subsystems. The numerical technique developed is then demonstrated using two exemplar applications: tissue engineering; and avascular tumour development. This allows verification that the technique is appropriate for many features of multiphase viscous flow modelling. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Multiphase viscous flow is a mathematical model that has been used in many areas of biology to model systems where two or more distinct phases co-exist in a given region. Each of these phases are assumed to behave like a viscous fluid, with problem specific inter-phase and intra-phase tractions specified by the application of interest. One example application of multiphase flow modelling is tissue engineering. In this application we model tissue that is grown externally from the subject by seeding a small number of healthy cells on a scaffold within a bioreactor, where nutrients for cell proliferation are provided via a culture medium flow. In this case the phases of the multiphase model are typically: (i) the cells and extra-cellular matrix; (ii) the culture medium; and (iii) the scaffold [1–4]. A second example is tumour development. For avascular tumour development [5,6] the phases are tumour cells and extra-cellular fluid. For vascular tumour development [7] the avascular tumour model may be adapted to include a blood vessel phase. Further examples of multiphase viscous flow models are: cell motility [8]; bloodflow [9]; and cartilage compression [10]. The system of Partial Differential Equations (PDEs) governing the above models consist of several components, and are discussed in detail in [1]. Each phase is modelled as a viscous fluid where the inertial term may be neglected. The volume fraction of each phase is governed by a non-linear hyperbolic PDE. This volume fraction is ⇑ Corresponding author. Address: Computing Laboratory, University of Oxford, Wolfson Building, Parks Road, Oxford OX1 3QD, United Kingdom. Tel.: +44 1865 273858; fax: +44 1865 273839. E-mail addresses: [email protected] (J.M. Osborne), Jonathan. [email protected] (J.P. Whiteley). 0045-7825/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2010.07.011

often dependent on various nutrient and messenger concentrations: these concentrations are modelled by elliptic PDEs. Further complexity is introduced in some cases, e.g. tumour development, by a time-dependent domain on which the solution must be calculated. In previous work [1–3,5–7,9,10], analytic techniques have been used to simplify the governing system of equations. Although this work has contributed significant insight into the underlying physical processes, many questions remain unanswered due to the simplifications that are necessary to apply these analytic techniques. In this study we present a technique for the numerical solution of the PDEs governing multiphase viscous flow. Our approach allows a more general coupling between phases than previous work [11]. This allows us to include more detail in our models, for example: (i) simultaneously taking account of several phases; (ii) allowing more general interactions between phases; and (iii) not requiring reduction to a one-dimensional model. The technique presented is based on a variety of finite element methods and the method of characteristics. We begin by writing down the governing PDEs in Section 2, and then deriving the numerical methods in Section 3. We then demonstrate the application of the technique to tissue engineering applications in Section 4 and avascular tumour development applications in Section 5. Finally, we discuss the merits of the scheme in Section 6. 2. The governing equations In this section we write down the governing equations for a multiphase viscous flow model, comprising N distinct phases and M distinct nutrient or messenger concentrations. These equations have been derived using the principles of conservation of mass

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

and momentum of each fluid phase, and conservation of mass of each nutrient or messenger concentration. See, for example, [1] for a discussion on the derivation of these equations. We will assume that the equations are valid on a time-dependent domain denoted by X(t), where T0 < t < T1. 2.1. Conservation of mass of each fluid phase 2.1.1. The governing PDEs Denoting the volume fraction, velocity and density of phase k by hk, uk and qk, conservation of mass of implies that hk satisfies the hyperbolic PDE

qk

  @hk þ $  ðhk uk Þ ¼ Sk ; @t

k ¼ 1; 2; . . . ; N;

x 2 XðtÞ;

ð1Þ

where Sk is a function that specifies net transfer of mass to phase k from other phases, and is specified by the model. These functions may depend on many different variables in the model and so we write, for k = 1, 2, . . ., N,

Sk ¼ Sk ðh1 ; . . . ; hN ; p1 ; . . . ; pN ; u1 ; . . . ; uN ; q1 ; . . . ; qN ; C 1 ; . . . ; C M Þ; where pk is the pressure within phase k, and Cl, l = 1, 2, . . ., M, is the concentration of nutrient or messenger l. We note that provided matter is only transferred between the phases, and not created or destroyed, we may write N X

Sk ¼ 0:

ð2Þ

3403

where pkl is the inter-phase pressure exerted on phase k by phase l (satisfying pkl = plk), and bkl is the viscous drag coefficient between phase k and phase l (satisfying bkl = blk, bkk = 0). The first term on the right hand side of Eq. (8) represents the inter-facial force exerted on phase k by phase l, with the second term representing the corresponding force exerted on phase l by phase k. The final term represents viscous drag between the phases. We now require model specific functional forms for rk, k = 1, 2, . . ., N, and pkl, k, l = 1, 2, . . ., N. Inter-phase pressures are given by

pkl ¼ p þ wkl ðh1 ; h2 ; . . . ; hN Þ;

k; l ¼ 1; 2; . . . ; N;

ð9Þ

where p is the overall mixture pressure and wkl is the inter-phase traction between phase k and phase l, and is specified by the application. Before specifying the stress tensor rk we first define the pressure of phase k, pk, by

pk ¼ p þ Rk ðh1 ; h2 ; . . . ; hN Þ þ

X

hl wkl ðh1 ; h2 ; . . . ; hN Þ;

k ¼ 1; 2; . . . ; N;

l–k

ð10Þ where Rk is a function prescribed by the application representing intra-phase pressure. Finally, using the assumption that each phase is a viscous fluid with dynamic shear viscosity lk and bulk viscosity kk, we may define the stress tensor as





rk ¼ pk I þ lk $uk þ ð$uk Þ> þ kk ð$  uk ÞI;

ð11Þ

where the gradient of a vector is defined by

k¼1

Making the assumption that there are no voids in the continuum allows us to deduce that the volume fractions sum to unity, allowing us to write N X

hk ¼ 1:

ð3Þ

k¼1

2.1.2. Boundary conditions and initial conditions Suitable boundary conditions for Eq. (1) are to specify hk, k = 1, 2, . . ., N, on all parts of the boundary for which uk flows into the computational domain. This may be written

hk ðxÞ ¼ hBk ðx; tÞ;

k ¼ 1; 2; . . . ; N;

x 2 @ Xk ðtÞ;

ð4Þ

hBk ;

for a prescribed function k ¼ 1; 2; . . . ; N, where oXk(t) is the inflow boundary for phase k satisfying uk  n < 0, where n is the outward-pointing unit normal to oX(t). Suitable initial conditions are to specify hk, k = 1, 2, . . ., N, at t = T0:

hk ðx; T 0 Þ ¼ h0k ðxÞ;

k ¼ 1; 2; . . . ; N;

x 2 XðT 0 Þ;

ð5Þ

where h0k ; k ¼ 1; 2; . . . ; N, is a prescribed function.

2.2.1. The governing PDEs Neglecting inertial effects, conservation of momentum of phase k may be written

 k

$  hk r þ

X

kl

f ¼ 0;

k ¼ 1; 2; . . . ; N;

ð6Þ

l–k

where rk is the stress tensor of phase k, and the divergence of a tensor r is defined by

ð$  rÞi ¼

X @ rij : @xj j

kl

ð12Þ

and I is the identity tensor. In common with other authors [1] we use kk = 2lk/3. 2.2.2. Boundary conditions Suitable boundary conditions for use in conjunction with Eq. (6) are specification of either the applied force or velocity in each spatial direction for each phase. Mathematically, this requires that we partition the boundary oX into two non-overlapping components for each phase: one component oXDk where Dirichlet boundary conditions are applied for phase k; and a second component oXNk where Neumann boundary conditions are applied for phase k. We may then write the boundary conditions as, for k = 1, 2, . . ., N:

uk ¼ u0k ;

for x 2 @ XDk ;

hk rk n ¼ sk ;

ð13Þ

for x 2 @ XNk ;

ð14Þ

where u0k and sk are prescribed functions.

The concentration of nutrient or messenger l is denoted by Cl, l = 1, 2, . . ., M. We assume that Cl satisfies the quasi-static diffusion equation

$2 C l ¼ Q l ðh1 ; h2 ; . . . ; hN ; C l Þ;

l ¼ 1; 2; . . . ; M;

ð15Þ

where Ql represents the net production of nutrient or messenger l. Boundary conditions are specification of either Cl or the flux of Cl at each point on the boundary. 2.4. The computational domain

ð7Þ

The term fkl denotes the force exerted on phase k by phase l, for k, l = 1,2,. . .,N, and is given by

f ¼ pkl hl $hk  plk hk $hl þ bkl hk hl ðul  uk Þ;

@ui ; @xj

2.3. Conservation of mass of each nutrient

2.2. Conservation of momentum of each fluid phase



ð$uÞij ¼

ð8Þ

At the beginning of Section 2 we explained that we allow the computational domain to vary with time. This is clearly necessary in some application areas, for example when modelling a growing tumour or crawling cell. We allow the domain to evolve with time by specifying the velocity of the boundary of the domain:

3404

dx ¼ v; dt

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

x 2 @ XðtÞ;

ð16Þ

where v is specified on each point of the boundary, possibly in terms of other dependent variables of the model. We note that we may allow one or more phase of the model to flow across the boundary oX(t), either into or out of the domain X(t). As such, the volume of the computational domain may vary with time without violating the condition that mass is not created or destroyed given by Eq. (2). We discuss this further when discussing application specific details for the examples given later in this paper.

3. The numerical scheme In this section we describe the numerical methods that we use to solve each of the equations introduced in Section 2. The relationships given by Eqs. (9)–(11) allow us to write all terms in the governing equations in terms of the variables u1, u2, . . ., uN, p, h1, h2, . . ., hN, C1, C2, . . ., CM. In order to compute each of these quantities numerically we uncouple the whole system of PDEs into the following subsystems:  N  1 equations of the form of Eq. (1) to compute h1, h2, . . ., hN1;  Eq. (3) to compute hN;  the velocity system (defined in Section 3.3) which is solved for u1, u2, . . ., uN and p;  M equations of the form of Eq. (15) to compute C1, C2, . . ., CM; and  updating the computational domain using Eq. (16), if required. These systems are then solved sequentially on each timestep until the final time is reached. Uncoupling in this way allows us to update each of the variables on each timestep in a systematic way: first we update the variables h1, h2, . . ., hN as these variables satisfy equations containing a time derivative. Using updated values of these variables then allows us to update the quantities u1, u2, . . ., uN, p, C1, C2, . . ., CM, and the boundary of the computational domain, based on the updated values of h1, h2, . . ., hN. We now describe how to solve each of the individual subsystems numerically. 3.1. The computational mesh

dhk 1 ¼ ðSk  hk $  uk Þ; dt qk dxk on ¼ uk : dt

ð17Þ ð18Þ

Writing each PDE as a system of ODEs at each point allows us to use the method of characteristics [12]: we may integrate Eq. (17) numerically to update hk, k = 1,2,. . .,N  1, at each point in space, whilst integrating Eq. (18) to update the point in space. In order to facilitate the solution of Eq. (17) we may introduce a small amount of artificial diffusion [13], this turns the first order hyperbolic equation to a second order parabolic equation of the form

dhk 1 ¼ ðSk  hk $  uk Þ þ Dk $2 hk ; dt qk

ð19Þ

where Dk is the artificial diffusion coefficient, which is taken to be zero unless otherwise indicated. Eqs. (18) and (19) can be solved by using the finite element technique with a linear or quadratic approximation to the solution on each element. See, for example, [15] for more details. The inclusion of such artificial diffusion has a negligible effect on the solution [13] but requires additional boundary conditions on the cell volume fraction; we impose @hk/@n = 0 on oX(t)n oXk (t), k = 1, 2, . . ., N, (i.e. where there is no inflow condition defined), where @/@n denotes the derivative normal to oX(t). When using the methods described above to solve Eq. (17) or (19), with (18) it is important to note that the presence of Eq. (18) will result in the functions hk, k = 1, 2, . . ., N  1, being computed at different points in space than at the previous timestep. On each timestep some extra work is required so that the updated solution is specified at the nodes of the finite element mesh. This may be achieved by interpolation on the finite element mesh. This will ensure that on each timestep an updated value of hk, k = 1, 2, . . ., N  1, is specified at each node of the finite element mesh. Finally, we use Eq. (3) to compute hN at the nodes of the mesh. 3.3. The velocity system Before describing how the velocity system is solved, we first list some mathematical definitions that underpin this component of the solution technique.

Before the variables are evolved on each timestep we assume that all dependent variables — u1, u2, . . ., uN, p, h1, h2, . . ., hN, C1, C2, . . ., CM — are known at the nodes of a finite element mesh. In the simulations performed later, we work in two spatial dimensions, and use triangular elements. The method presented here may, however, be adapted for use with quadrilateral meshed in two dimensions, or both tetrahedral or hexahedral elements in three dimensions. As discussed above, we use a finite element mesh that partitions the computational domain into non-overlapping triangular elements. We assume two triangles intersect along a whole edge, at a single node, or not at all: i.e. there are no hanging nodes in the mesh. If we know the value of a quantity at each node of the mesh we may approximate the value of that quantity at any point within the mesh: we simply identify which element contains the point of interest, and then interpolate the nodal values of that element to compute the approximation to that value.

The following subsets of H1(X) for scalar valued functions are then defined by

3.2. Updating the volume fraction of each phase

Finally, the compound subspaces H1Ek and H1Ek0 for vector valued functions are defined by

For k = 1, 2, . . ., N  1 we may re-write Eq. (1) as the following system of ordinary differential equations (ODEs) at each point x 2 X(t):

n H1Ek ¼ uk 2 H1 ðXÞ  H1 ðXÞ  H1 ðXÞ : uk satisfies o specified Dirichlet boundary conditions ;

3.3.1. Function spaces and other preliminaries The spaces L2(X) and H1(X) are defined by

 Z u: u2 dV < 1 ; X   @u 1 2 L2 ðXÞ; i ¼ 1; 2; 3 : H ðXÞ ¼ u 2 L2 ðXÞ : @xi

L2 ðXÞ ¼



n H1E ðXÞ ¼ u 2 H1 ðXÞ : u satisfies specified Dirichlet o boundary conditions ; n H1E0 ðXÞ ¼ u 2 H1 ðXÞ : u satisfies zero Dirichlet o boundary conditions :

3405

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

n



sk ¼ lk $uk þ ð$uk Þ

H1Ek0 ¼ uk 2 H1 ðXÞ  H1 ðXÞ  H1 ðXÞ : uk

o satisfies zero Dirichlet boundary conditions :

r:s¼

rij sij :

ð20Þ

3.3.2. The weak form of the velocity system In D spatial dimensions, Eq. (6) provides us with a system of PDEs of size ND. However, we wish to solve these equations for u1, u2, . . ., uN, p: a total of ND + 1 quantities. Noting that, in Section 3.2, we only used N  1 of the N equations that comprise Eq. (1), we may still use one equation derived from this system. We will use the sum of this equation over all phases, and will now show that this sum may be considered analogous to the incompressibility condition when solving the Stokes equations governing flow of a single viscous fluid with negligible inertial effects [14]. Summing Eq. (1) over all phases yields

qk

k¼1

  X N @hk Sk : þ $  ðhk uk Þ ¼ @t k¼1

Application of: (i) the no voids condition given by Eq. (3); (ii) the condition that there is no net creation of mass given by Eq. (2); and (iii) assuming that qk is identical for each phase allows us to write

$

N X

! ¼ 0:

hk uk

ð21Þ

k¼1

Defining the mixture velocity U by

U¼

N X

hk uk ;

k¼1

implies that Eq. (21) reduces to $  U ¼ 0, which is of the same form as the incompressibility constraint associated with the Stokes equations. We may now augment Eq. (21) to the system given by Eq. (6) so that there are now ND + 1 equations for the ND + 1 quantities u1, u2, . . ., uN, p. We now turn our attention to deriving the weak form of the system described by Eqs. (6) and (21). For k = 1, 2, . . ., N we take the scalar product of Eq. (6) with wk 2 H1Ek0 , and integrate over X. We have chosen to use wk 2 H1Ek0 so that wk = 0 on oXDk. Using this observation, application of the divergence theorem, and of the Neumann boundary condition (Eq. (14)) yields, after some manipulation:

Ak ðu1 ; u2 ; . . . ; uN ; wk Þ  Bk ðp; wk Þ ¼ C k ðwk Þ;

8wk 2 H1Ek0 ðXÞ; ð22Þ

where

Ak ðu1 ; u2 ; . . . ; uN ; wk Þ ¼

Z

hk ð$wk Þ : sk 

X

Bk ðp; wk Þ ¼ Z C k ðwk Þ ¼

Z X

X

Z X

N X

Bk ðp; wk Þ ¼

k¼1

N X

C k ðwk Þ:

ð27Þ

k¼1

The weak formulation of Eq. (21) is derived by multiplying this equation by q 2 L2(X) and integrating over X. This yields N X

Bl ðq; ul Þ ¼ 0;

8q 2 L2 ðXÞ:

ð28Þ

l¼1

Writing

Aðu1 ; u2 ; . . . ; uN ; w1 ; w2 ; . . . ; wN Þ ¼

N X

Ak ðu1 ; u2 ; . . . ; uN ; wk Þ; ð29Þ

k¼1

Bðp; w1 ; w2 ; . . . ; wN Þ ¼

N X

Bk ðp; wk Þ;

ð30Þ

k¼1

Cðw1 ; w2 ; . . . ; wN Þ ¼

N X

C k ðwk Þ;

ð31Þ

k¼1

we see that the weak formulation of Eqs. (6) and (21) may be written: find u1 2 H1E1 ðXÞ; u2 2 H1E2 ðXÞ; . . . ; uN 2 H1EN ðXÞ; p 2 L2 ðXÞ, such that, for all wk 2 H1Ek0 ðXÞ; k ¼ 1; 2; . . . ; N; q 2 L2 ðXÞ

Aðu1 ; u2 ; . . . ; uN ; w1 ; w2 ; . . . ; wN Þ  Bðp; w1 ; w2 ; . . . ; wN Þ ¼ Cðw1 ; w2 ; . . . ; wN Þ;

ð32Þ

Bðq; u1 ; u2 ; . . . ; uN Þ ¼ 0:

ð33Þ

For more details on the derivation of this weak formulation see Appendix A. We note the similarity of Eqs. (32) and (33) with the weak formulation of the Stokes equations [15]: this will be used in the next section to guide our choice of finite element method. 3.3.3. The finite element solution of the velocity system Given the similarity of Eqs. (32) and (33) with the weak formulation of the Stokes equations, we will use an established finite element method that has been used for the Stokes equations, namely a quadratic approximation on each element to each component of each phase velocity uk, k = 1, 2, . . ., N, and a linear approximation on each element to the pressure p. For details on the application of this finite element method see Elman et al. [15]. We denote by V hEk ðXÞ and V hEk0 ðXÞ the finite dimensional subspaces of H1Ek ðXÞ and H1Ek0 ðXÞ that consist of quadratic functions on each element that are continuous across element boundaries. We also denote by Wh the finite dimensional subspace of L2(X) that consists of linear functions on each element that are continuous across element boundaries. Denoting the finite element solution by uhk ; k ¼ 1; 2; . . . ; N, and ph, this solution satisfies: find uh1 2 V 1E1 ðXÞ; uh2 2 V 1E2 ðXÞ; . . . ; uhN 2 V 1EN ðXÞ; ph 2 W h , such that, for all wk 2 V 1Ek0 ðXÞ; k ¼ 1; 2; . . . ; N; q 2 W h ,

ð34Þ

ð23Þ

Bðq; uh1 ; uh2 ; . . . ; uhN Þ ¼ 0:

ð35Þ

ð24Þ

This finite element method is used to compute the numerical approximation to the phase velocities and pressure on each timestep.

l–k

wkl ðhl $hk  hk $hl Þ  wk dV;

X l–k k

Ak ðu1 ; u2 ; . . . ; uN ; wk Þ 

¼ Cðw1 ; w2 ; . . . ; wN Þ;

l–k

X

ð26Þ

Aðuh1 ; uh2 ; . . . ; uhN ; w1 ; w2 ; . . . ; wN Þ  Bðph ; w1 ; w2 ; . . . ; wN Þ bkl hk hl ðul  uk Þ  wk dV;

p$  ðhk wk ÞdV; Z X wk  sk dS þ hk ðRk þ hl wkl Þ$  wk dVþ

@ XNk

N X k¼1

i;j

N X

þ kk ð$  uk ÞI :

Summing Eq. (22) over all phases k = 1, 2, . . ., N we obtain, for all wk 2 H1Ek0 ðXÞ; k ¼ 1; 2; . . . ; N,

We will also require the tensor product of two tensors r and s, written r:s, defined by

X

>

and s is the viscous component of the stress tensor given by

ð25Þ

3.4. Solution of Eq. (15) As h1, h2, . . ., hN have all been computed using the methods described in Section 3.2, each equation given by Eq. (15) is an elliptic

3406

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

PDE for the single dependent variable Cj, j = 1, 2, . . ., M. These equations may be solved using the finite element technique with a linear approximation to the solution on each element. See, for example, [15] for more details on the application of this technique. 4. Tissue engineering applications The first exemplar application we consider is that of a multiphase viscous flow model of tissue growth, external to the body, within a perfusion bioreactor. We begin by giving a brief overview of this application in Section 4.1 before writing down the model specific functions required to fully pose the model in Section 4.2. We then validate our software in Section 4.3 by verifying that our solutions converge to a simpler model in an appropriate limit, and then demonstrate that the numerical technique may be used for more general simulations in Section 4.4. Finally, in Section 4.5, we discuss how a quantity of interest varies with the spatial mesh and timestep used. 4.1. The bioreactor We consider the mathematical model of the compression perfusion bioreactor geometry described by El Haj and co-workers [16– 18]. This bioreactor provides both perfusion and macroscale compression to a culture of human-derived bone cells contained within a biodegradable poly L-lactic acid (PLLA) scaffold. The system comprises a porous scaffold approximately 8 mm in diameter and 4 mm in height, with a mean porosity of 97%, seeded with osteoblast like cells, which is placed within a culture medium filled cylinder. A flow is driven across the scaffold (perfusion at rate 0.1 ml/ min due to a peristaltic pump), allowing the mechanical stimulation of cells via pressure and/or fluid shear. In addition there is a piston in the bioreactor which can compress the scaffold with strains of up to 1.5% at a rate of 1 Hz. We do not consider compression component of the bioreactor here, however the mathematical model given here could be extended to include the effects of the piston. This system is illustrated schematically in Fig. 1. 4.2. Model equations We now write down the model specific functions and boundary conditions that are required to close the model. Unless otherwise stated in our examples we work with non-dimensionalised equations and variables, hence no units are usually given. Where appropriate we translate the presented results back into dimensional form, to enable comparisons to be made with experimental data thus demonstrating that we are using realistic parameter values. 4.2.1. The phases in the model The tissue engineering model that we use contains: (i) a cell phase, comprising cells and extra-cellular matrix; (ii) a fluid phase, comprising culture medium and extra-cellular fluid; and (iii) a solid scaffold phase that cannot be modelled as a viscous fluid. Variables associated with the scaffold phase are denoted with the subscript s. We assume that: (i) the scaffold volume fraction

Fig. 1. The bioreactor system of [17]. (Figure adapted from [3].)

is a specified time-independent function hs(x); and (ii) the scaffold is fixed and so us = 0. We therefore do not need to solve any equations to calculate these quantities. The other two phases are modelled as viscous fluids, and the numerical techniques described in Section 3 may be used to compute the variables associated with these phases. We note, however, that we still allow the effect of intra-phase pressures exerted by the solid scaffold phase to be included in the model. This model is described in more detail in O’Dea et al. [3] and Osborne et al. [19]. To use Eqs. (9) and (10) we need to specify intra-phase and inter-phase pressure functions. We use those used in previous work [1,3,19]. We use a subscript n to represent the cell phase, and a subscript w to represent the fluid phase. The intra-phase and inter-phase functions are all zero, with the exception of Rn and wns, which are defined by



Rn ¼ hn m þ

 da hm n ; hm w

and wns ¼ v þ

db hm n ; hm w

ð36Þ

for constants m, v, da, db, m > 0 whose physical interpretation and values are given in Table 1. Some of these parameters are scaled using another parameter h which is the width of the bioreactor. This scaling allows us to take the limit h ? 0, in which case a one-dimensional limit exists [3]: this limit is used to verify the accuracy of our software. 4.2.2. Boundary and initial conditions The bioreactor described in Section 4.1 is represented by a fixed two-dimensional channel X: 0 6 x 6 1, 0 6 y 6 h where the ends x = 0, 1 are open, and the sides y = 0, h are fixed channel walls. Appropriate conditions on the channel walls (y = 0, h) are no-slip and no-penetration (of cells or culture medium). To model perfusion, a pressure-driven flow is imposed via up- and downstream pressures, PU, and PD. This leads to normal stress conditions on each phase at x = 0, 1 in which we assume that each phase bears its share of the stress in proportion to its volume fraction. We further assume that the flow both upstream and downstream is fully developed. These boundary conditions are the same as used in [19] and are illustrated in Fig. 2, in which n is the outward-pointing unit normal to the boundary oX of X. As discussed above, we assume that the scaffold volume fraction is constant in time. In our simulations the scaffold volume fraction, hs is defined by

 hs ðx; y; tÞ ¼

0:03 for a 6 x 6 b; 0

otherwise;

ð37Þ

where 0 < a < b < 1 are the edges of the scaffold. The values of these parameters used in this work are a = 0.25 and b = 0.75. Along with the scaffold volume fraction, initial conditions for the cell volume

Table 1 Table of general parameters in the bioreactor model. Parameters that are fixed in the simulations PD Downstream applied pressure PU Upstream applied pressure D Artificial diffusion coefficient v Scaffold affinity parameter m Cell aggregation parameter da Intra-phase repulsion parameter db Inter-phase repulsion parameter m Repulsion scaling parameter ln Ratio of dimensional viscosities ln/lw b Inter-phase drag coefficient a Left boundary of scaffold b Right boundary of scaffold

0.1/h2 0.3/h2 0.001 0.05/h2 0.05/h2 0.05/h2 0.05/h2 1 1.3 10 0.25 0.75

Parameters that are varied in the simulations h Width of the reactor (aspect ratio)

h 2 (0, 1]

3407

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417 Table 2 Table of parameters for cell proliferation function in the bioreactor model.

Fig. 2. The two-dimensional domain of the bioreactor, X, and associated boundary conditions.

fraction, hn0 ðx; yÞ ¼ hn ðx; y; t ¼ 0Þ, need to be defined. As hn + hw + hs = 1 initial conditions for hw are not required. Suitable initial conditions for the simulations carried out here are given in Sections 4.3 and 4.4. In these simulations the computational domain is fixed, and so we set v = 0 in Eq. (16). 4.2.3. The cell proliferation functions Finally we close the model by defining the proliferation functions, Sn and Sw, required to solve Eq. (1). Following [3,19], we demonstrate the application of our numerical algorithm for four mechanisms for cell proliferation and growth:

Constant proliferation Sn(hn) km Rate of mitosis kd Rate of apoptosis

1.0 0.1

Cell density-dependent proliferation Sn(hn) hn1 Density threshold Density threshold hn2 Rate of mitosis k1h Rate of mitosis and ECM deposition k2h Rate of apoptosis kdh

0.4 0.6 0.8 1.0 0.1

Cell pressure-dependent proliferation Sn(hn, pn) pn1 Pressure threshold pn2 Pressure threshold Rate of mitosis k1p

0.15/h2 0.22/h2 0.8

k2p

Rate of mitosis and ECM deposition

1.0

kdp

Rate of apoptosis

0.1

Fluid phase shear stress dependent proliferation Sn(hn, sw) Stress threshold Stress threshold Rate of mitosis k1s Rate of mitosis and ECM deposition k2s Rate of apoptosis kds

sw1 sw2

0.025h 0.125h 0.5 1.0 0.1

 uniform proliferation, with Sn = Sw = (km  kd)hn;  proliferation dependent on the cell density (representing contact inhibition-regulated growth), where Sn = Sw = j(hn)hn;  proliferation dependent on the cell density and the cell pressure, where Sn = Sw = j(pn)hn; and  proliferation dependent on the cell density and the flowinduced shear stress, sw, given by [14]:

sw ¼

 2  2 !12 @uw @v w þ ; @y @x

ð38Þ

where Sn = Sw = j(sw)hn. For the last three cases above we define a net tissue growth rate j(/) such that Sn = j(/)hn where / = hn, pn, sw denotes density-, pressureand culture medium shear stress-mediated growth. For each of these regulatory stimuli, three distinct stages in the cell population’s behaviour have been identified: a quiescent phase (proliferation and ECM deposition are reduced), a proliferative phase (proliferation and ECM deposition are elevated) and an apoptotic phase. In each of these stages j(/) takes a different value. The thresholds separating these stages are denoted /1 and /2, where 0 < /1 < /2. A physical description of the parameters appearing in the relationships above, and the values used in this study, are given in Table 2. For density-dependent growth and pressure-dependent growth, a step function representation is employed such that the net growth rate takes the values k1/ for 0 6 / < /1, k2/ for /1 6 / 6 /2, and kd for / > /2. In contrast, to obtain a stable numerical solution, in the case of culture medium shear stress-mediated growth, the discontinuities in j(sw) at sw1 and sw2 must be smoothed (see [3,19] for more details). The functional forms employed for j in the above regimes are specified below:









jð/Þ ¼ k1/ Hð/1  /Þ þ k2/ Hð/  /1 Þ  k2/ þ kd Hð/  /2 Þ ; / ¼ hn ; p n ;     k1s  k1s  tanh g s  sw1  1 2     k2 þ kds  tanh g sw  sw2  1  kds ;  s 2

ð39Þ

jðsw Þ ¼ 

ð40Þ

in which H(/) is the Heaviside function. For simulations presented in this section a value of g = 50 is used in Eq. (40). Note that in

Fig. 3. An example mesh used for the tissue engineering simulations.

the limit g ? 1 we recover a shear-dependent step function response. 4.2.4. Numerical parameters used All simulations presented in this section are produced using the following numerical parameters. When running simulations on the domain X = [0, 1]  [0, 1] we use a mesh consisting of 30 elements in each direction. This results in 31  31 = 961 nodes for linear basis functions (used to represent pressure), and 61  61 = 3721 nodes for the quadratic basis functions (used to represent volume fractions and phase velocities). For meshes on the domain X = [0, 1]  [0, h], the number of elements in the y direction is adjusted: for X = [0, 1]  [0, 0.01] and X = [0, 1]  [0, 0.1] we use 40 elements in the x direction and 5 elements in the y direction; and for X = [0, 1]  [0, 0.2] we use 40 elements in the x direction and 10 elements in the y direction. The nodes are concentrated near the scaffold (a 6 x 6 b) to capture the variables as efficiently as possible as variation is low away from the scaffold. An example mesh for the domain X = [0, 1]  [0, 0.1] is shown in Fig. 3. All simulations use a timestep of 1  104. In order to stabilise the numerical solutions artificial diffusion of value Dn = Dw = D = 0.001 is used. 4.3. Example results: convergence of solutions to a 1D model The thin-film approximation of the model is a simplification of the bioreactor model, and is found by considering the bioreactor to

3408

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

be long and thin, i.e. the limit h ? 0. This allows analytic simplification of the two-dimensional problem to a one-dimensional problem [3] for which numerical packages may be used to compute the numerical solution. Note that this approximation is referred to as the long-wavelength limit in [3]. In this section we present results of simulations of bioreactors with varying aspect ratios in order to ascertain convergence of the solutions to the thin-film approximation, thus giving some verification of the correctness of our software. Following [19], we consider the evolution of the tissue construct composition from the following initial state:

hn0 ðx; yÞ ¼ 0:25½tanhð75ðx  0:4ÞÞ  tanhð75ðx  0:6ÞÞ:

ð41Þ

This condition, along with Eq. (37), represents a small population of cells (of approximate width 0.2) distributed within a scaffold (of width b  a) of uniform porosity 97%. Fig. 4 shows simulations with uniform cell proliferation. Cell volume fractions at time t = 1 are presented for bioreactors of width h = 0.2 in Fig. 4(a), h = 0.1 in Fig. 4(b), and h = 0.01 in Fig. 4(c). In Fig. 4(d) we plot the cell volume fraction averaged over the width of the bioreactor  hn ðx; tÞ, given by

hn ðx; tÞ ¼

Z

h

!, hn ðx; y; tÞdy

h;

ð42Þ

0

at time t = 1 for h = 0.2, 0.1, 0.01, together with the value of  hn computed from the thin-film approximation. We see that our numerical results converge to the thin-film approximation of the model (shown as the dotted line) in the limit h ? 0. As such, we have confidence in our software.

4.4. Example results: comparison of mechanotransduction mechanisms In Section 4.3 we demonstrated convergence of the two dimensional solutions to the solutions of a simplified model. We now consider a more complex example. One common tissue engineering experiment is to seed cells onto the scaffold by placing the scaffold into a cell culture allowing the cells to penetrate the scaffold. This may be modelled by using an initial cell volume fraction of

!! pffiffiffi 4 1 2 21 pffiffiffi  tanh 50 x  hn0 ðx; yÞ ¼ 10 10 4 2 ! !! pffiffiffi 2 2þ1 pffiffiffi  tanh 50 x  4 2 !! pffiffiffi 21  tanh 50 y  pffiffiffi 2 2 !!! pffiffiffi 2þ1  tanh 50 y  pffiffiffi ; 2 2

ð43Þ

for 0.25 6 x 6 0.75, and hn0 ðx; yÞ ¼ 106 , otherwise. This initial condition is shown in Fig. 5. In order to illustrate the general applicability of our numerical scheme in Fig. 6 we present simulations using all mechanotransduction mechanisms discussed in Section 4.2.3. These results are shown in Fig. 6(a) for uniform proliferation, Fig. 6(b) for cell density-dependent growth, Fig. 6(c) for cell pressure-dependent growth, and Fig. 6(d) for fluid phase shear stress dependent growth. In all figures we plot the solution at t = 0.5. It is clear from Fig. 6 that our solution techniques (and software) can handle the

Fig. 4. Simulations with uniform cell proliferation. Cell volume fractions at time t = 1 for bioreactors of varying widths: (a) h = 0.2; (b) h = 0.1; and (c) h = 0.01. (d) Comparison of averages across bioreactors of widths: h = 0.2, ‘‘--”; h = 0.1, ‘‘- - -”; h = 0.01, ‘‘—”; and the thin-film approximation, ‘‘  ”. Parameter values are given in Tables 1 and 2.

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

3409

and  hn that arise as a consequence of different cell proliferation specifications. 4.5. Example results: dependence of the solution on spatial mesh and timestep

Fig. 5. Initial cell volume fraction, hn0 , for simulations of tissue growth in the bioreactor.

two dimensional phenomena that are modelled by these example simulations. We now consider some processing of the numerical solution that would be of interest to experimental workers in this field. The cell volume fraction averaged across the bioreactor is given by Eq. (42). The cell yield,  hn ðtÞ is given by

h ¼ n

Z

h

y¼0

Z

1

hn ðx; y; tÞdxdy:

ð44Þ

x¼0

In Fig. 7(a) we plot  hn as a function of x at t = 0.5, and in Fig. 7(b) we plot  hn as a function of time. In these plots, the solid line represents uniform cell proliferation, the broken line represents cell densitydependent growth, the dash-dot line represents cell pressuredependent growth and the dotted line represents fluid shear stress dependent growth. This figure demonstrates that our numerical technique may be used to demonstrate the differences in both  hn

In this section we investigate how one quantity of interest that may be computed from the numerical solution — namely the yield given by Eq. (44) — depends on the spatial mesh and timestep. In Fig. 8(a) we investigate how the yield varies with spatial mesh. We use the parameters given in Section 4.2.4, with the exception of the nodal spacing in the computational mesh. In these computations, we use a uniform mesh, with spacing dx between nodes in the x-direction, and spacing 0.1dx between nodes in the y-direction. We plot the yield for: dx = 0.1 (solid line); dx = 0.05 (broken line); dx = 0.025 (dot-dashed line); and dx = 0.0125 (dotted line). We see that there is very little difference between the computed yield for these simulations. In Fig. 8(b) we investigate how the yield varies with timestep. We use the mesh corresponding to dx = 0.025 in Fig. 8(a), and compute the yield for a timestep of: 2  104 (solid line); 1  104 (broken line); and 5  105 (dotted line). In common with the simulations shown in Fig. 8(a) we see that the timestep of 104 used in our simulations for these model parameters is appropriate. We did not encounter any problems with numerical stability when the spatial mesh and timestep were varied. We emphasise here that the convergence of the yield shown in Fig. 8 is only valid for the model parameters used in this study. Should different model parameters be used, it is then likely that the solution may vary on different lengthscales and different timescales. Under these circumstances a different spatial mesh and a different timestep may be appropriate.

Fig. 6. Cell volume fraction at time t = 0.5 for: (a) uniform proliferation, (b) cell density-dependent growth; (c) cell pressure dependent growth; and (d) fluid shear stress dependent growth. Simulations are in a bioreactor of width h = 1.0 for time t = 0–0.5, with parameters as in Tables 1 and 2 and initial conditions shown in Fig. 5.

3410

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

Fig. 7. Functions of cell volume fraction at time t = 0.5. (a) Shows the cell volume fractions averaged across the bioreactor and (b) shows the cell yield over time for: uniform proliferation, ‘‘—”; cell density-dependent growth, ‘‘- - -”; cell pressure-dependent growth, ‘‘--”; and fluid phase shear stress dependent growth, ‘‘  ”. Simulations are in a bioreactor of width h = 1.0, with parameters as in Tables 1 and 2 and initial conditions shown in Fig. 5.

Fig. 8. (a) Convergence of the yield with spatial mesh and (b) convergence of the yield with timestep.

5. Avascular tumour development applications Our second exemplar application is a multiphase model of avascular tumour development. The model that we use is a generalisation of the one presented in [20]. The main properties of the model are:  the tumour consists of two phases, cells and extra-cellular fluid;  cell proliferation is dependent on nutrients which diffuse throughout the tumour; and

 the model is posed on a growing domain, where the movement of the boundary of the domain, and thus the evolution of the whole domain, is a consequence of the variables in the tumour model and not known a priori. 5.1. Model equations We again represent variables associated with the cell phase with a subscript n, and variables associate with the fluid phase with a subscript w. For this model, all the intra-phase and

Table 3 Parameters used in the tumour growth model. Parameters for nutrient equation Q0 Nutrient up-take parameter

0.5

Parameters for volume fraction equation s1 Cell proliferation s2 Cell proliferation s3 Cell proliferation s3 Cell proliferation

10 0.5 0.5 10

parameter parameter parameter parameter

Parameters for velocity system Cell attraction threshold hmin n

hn r p

ln lw b

Natural cell density Inter-phase scaling parameter Inter-phase scaling parameter Viscosity of cell phase Viscosity of fluid phase Inter-phase drag coefficient

0.8 0.8 1 1 1 0.1 1

Fig. 9. The initial mesh used for simulations of growing tumours.

3411

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

inter-phase pressures required by Eqs. (9) and (10) are assumed to be zero with the exception of Rn, which is given by [5,20]

Rn ¼

8 < 0; :c

jhn hn jr1 q hw



for 0 6 hn < hmin n ;   min hn  hn ; for hn 6 hn < 1; hmin n

hn

ð45Þ

for positive constants q, r, c and 0 < 6 < 1. The parameters q; r; c; hmin and hn depend on the characteristics of the cell types n being studied. The functions representing conversion of one phase to another are given by [5,20]

Sn ¼ Sw ¼

s0 C s2 þ s3 C hn hw  hn ; a þ s1 C 1 þ s4 C

where C represents nutrient concentration. Finally, nutrient consumption is given by

Q C ¼ Q 0 Chn : The parameter values used in all of these functions are taken from [5] and are listed in Table 3. 5.1.1. Boundary and initial conditions We now prescribe initial and boundary conditions for the tumour model. Initially the tumour has boundary oX0, and cell volume fraction hn0 ðxÞ, and so

@ Xð0Þ ¼ @ X0 ;

and hn ðx; 0Þ ¼ hn0 ðxÞ;

for x 2 @ Xð0Þ:

ð46Þ

Fig. 10. Tumour growth with uniform applied nutrient. Cell volume fraction, hn, throughout the tumour for times: (a) t = 0; (b) t = 25; (c) t = 50; (d) t = 75; and (e) t = 100, corresponding to 14 days. (f) The size of the tumour over time, in the x dimension, ‘‘—”, and in the y dimension, ‘‘- - -”. All parameters are as given in Table 3.

3412

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

The boundary of the tumour, oX(t), evolves over time. We assume that the tumour boundary moves at the velocity of the cell phase, and so we may write

dx ¼ un ; dt

for x 2 @ XðtÞ;

8t:

ð47Þ

Specifying the evolution of the boundary of the domain in this way allows the volume of the domain to change with time. In this application, which models a growing tumour, the velocity of the cell phase will be such that the volume of the tumour grows with time due to Eq. (47). To maintain consistency with Eq. (2), the equation that models mass neither being destroyed or created, this extra volume is created by fluid phase flowing into the domain with fractional concentration 1  hn. In previous work [5,20], the nutrient concentration outside the tumour has been assumed to be constant. Our approach allows a more general function that can vary in time and space, CX(x, t) to

be used. The boundary condition for the nutrient equation is therefore

Cðx; tÞ ¼ C X ðx; tÞ;

for x 2 @ XðtÞ;

8t:

ð48Þ

For this work we assume that the tumour is growing in free suspension (i.e. there is no external medium present), and assume zero stress on the boundary of the tumour. The boundary condition for the velocity system, is therefore given by

rn  n ¼ rw  n ¼ 0; for x 2 @ XðtÞ; 8t;

ð49Þ

where n is the outward-pointing unit normal to the boundary oX(t). 5.1.2. Numerical parameters used All simulations presented in this section are produced using the following numerical parameters. Simulations were run using an

Fig. 11. Tumour growth with uniform applied nutrient: (a) the cell volume fraction, hn; (b) the nutrient distribution, C; (c) and (d) the x and y component of the velocity of the cell phase, un and vn, respectively; (e) the radial component of the velocity of the cell phase, unr ; and (f) the bulk pressure, p, throughout the tumour. All plots are at time t = 75. All parameters are as given in Table 3.

3413

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

initial circular mesh consisting of 1594 elements which is shown in Fig. 9. We use a higher density of elements closer to the boundary of the tumour: this is to capture the more rapidly varying features of the solution in this region. The initial mesh is updated at each timestep by moving all nodes in the mesh by linear interpolation of the movement of the boundary. In addition the domain can be remeshed in order to maintain the quality of the mesh, for example should an element become larger than some threshold value, however, for these simulations no remeshing is undertaken. All simulations use a timestep of 5  102. For these simulations no artificial diffusion is added, and so therefore Dn = Dw = 0. 5.2. Example results: uniform growth We begin by simulating radially symmetric growth, allowing us to validate our software by comparing with the one-dimensional radial symmetric solution. This requires us to use the following non-dimensional boundary condition for C:

Cðx; tÞ ¼ 1;

for x 2 @ XðtÞ;

8t:

ð50Þ

This is uniform around the boundary of the tumour, and the problem is therefore radially symmetric. We show the solution in Figs. 10–12. Fig. 10(a)–(e) show the region occupied by the tumour at times t = 0, 25, 50, 75, 100, respectively (corresponding to growth over 14 days), where the degree of shading inside the tumour represents the cell volume fraction. Fig. 10(f) shows the radius of the tumour as a function of time. The resulting growth shown in Fig. 10 is clearly radially symmetric, as expected. As the tumour becomes large enough a necrotic core develops, where the cell volume fraction is reduced due to low concentration of nutrient, as shown by the lighter, interior regions in Fig. 10(d) and (e). In addition Fig. 10(f) shows that the spheroid grows from 50 lm to about 600 lm in diameter over a dimensional time of 14 days (as one unit of dimensionless length corresponds to 25 lm). This agrees exceptionally well with the experimental, free suspension, growth curves contained in [21] and also the simulated growth curves published in [5]. Note that the mathematical model contained in [5] is a simplification of the one presented here, where the tumour is assumed to be onedimensional and composed of viscous cells and inviscid extra-cellular fluid, therefore, direct quantitative comparison with these results is not possible. In Fig. 11 we present the variables of the same simulation at a representative time (t = 75). It is clear that all variables are radially symmetric. Fig. 11(e) demonstrates that the growth of the tumour is due to the proliferation of cells in an outer proliferating rim. On

this outer proliferating rim cells are flowing outwards ðunr > 0Þ and on the inside edge of the rim cells are flowing towards the centre of the tumour ðunr < 0Þ. In the central necrotic region there is no flow ðunr ¼ 0Þ. From Fig. 11(e) we see that the non-dimensional velocity of the rim of the tumour is approximately 0.1, therefore, the dimensional velocity is approximately 2  1010 m s1 corresponding to a growth rate of about 36 lm per day, which is the growth rate seen in Fig. 10(f) and [21]. We show plots of hn, C, un and p through the cross-section y = 0 at time t = 75 in Fig. 12(a), and plots of hn, C, vn and p through the cross-section x = 0 in Fig. 12(b). Noting that the first plot shows un whilst the second of these plots shows vn, we see that these plots are identical, demonstrating the radial symmetry of this problem. 5.3. Example results: non-homogeneous growth We now consider the growth of a tumour exposed to a nonhomogeneous external nutrient concentration, allowing us to demonstrate the application of our numerical technique to problems with marked two-dimensional features. The following non-dimensional boundary condition for C is applied:

Cðx; tÞ ¼

y þ 20 ; 20

for x 2 @ XðtÞ;

8t:

ð51Þ

This corresponds to a non-dimensional nutrient concentration of 2 at y = 20 decreasing linearly to zero at y = 20. Fig. 13(a)–(e) show the region occupied by the tumour at times t = 0, 25, 50, 75, 100, respectively, where the degree of shading inside the tumour represents the cell volume fraction. Fig. 13(f) shows the radius of the tumour as a function of time, where the solid line represents the radius in the x-direction, and the broken line represents the radius in the y-direction. Fig. 14 shows the solution of various dependent variables in the model at a non-dimensional time t = 75. We show plots of hn, C, un and p through the cross-section y = 0 at time t = 75 in Fig. 15(a), and plots of hn, C, vn and p through the cross-section x = 0 in Fig. 15(b). From Figs. 13–15 it is clear that the resulting growth is no longer radially symmetric, in fact the tumour ‘‘grows up” the nutrient gradient causing the tumour cells to proliferate towards the regions of higher nutrient concentration seen in Fig. 14(b). This growth is only due to the nutrient dependent growth mechanism affecting the velocity of the cell phase (Fig. 14(c)–(e)) and not a more complex signalling pathway. Fig. 13(f) shows the tumour size over time, the tumour grows unbounded in all directions except negative y, this is because the non-dimensional nutrient level at y = 5 is so small that the cells cease to proliferate.

Fig. 12. Cross sections of various quantities taken from Fig. 11. (a) Cross section along the line y = 0 (b) cross-section along the line x = 0 for a uniform external applied nutrient.

3414

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

Fig. 13. Tumour growth with non-uniform external applied nutrient. Cell volume fraction, hn, throughout the tumour for times: (a) t = 0; (b) t = 25; (c) t = 50; (d) t = 75; and (e) t = 100 corresponding to 14 days. (f) The size of the tumour over time, in the x dimension, ‘‘—”, and in the y dimension, ‘‘- - -”. All parameters are as given in Table 3.

6. Discussion In this study we have presented a numerical technique for computing the numerical solution of the equations governing multiphase viscous flow, and applied this technique successfully to exemplar applications from tissue engineering and avascular tumour development. These examples were chosen because they contain all the difficulties inherent with computing a numerical solution of these equations, namely: (i) the solution of different PDE types; and (ii) a time-dependent computational domain. The numerical solutions all converged when the spatial mesh was made finer and the timestep was made shorter, and no problems with numerical stability or dissipation were encountered. Based on these results, we see no reason why these techniques cannot be applied to other multiphase viscous flow models of the form considered in this study.

Future challenges for the development of the numerical framework include incorporating more physical detail, such as: implementation of poro- and visco-elastic multiphase models; modelling the degradation and compression of the scaffold phase in tissue engineering; and coupling the multiphase models to discrete models of vascular formation or interacting tissues. Efficient computation of an accurate numerical solution of the governing equations will allow the models being developed in several application domains to be validated by comparison with experimental data. In Section 4.5 we explained that it is not possible to prescribe explicitly the granularity of the spatial grid or the timestep that should be used for a general model parameter set. The development of the techniques presented here to allow adaptive spatial meshing and choice of timestep will allow computation of a numerical solution to within a user-specified tolerance, with an appropriate spatial mesh and timestep chosen automatically.

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

3415

Fig. 14. Tumour growth with non-uniform external applied nutrient: (a) the cell volume fraction, hn; (b) the nutrient distribution, C; (c) and (d) the x and y component of the velocity of the cell phase, un and vn, respectively; (e) the radial component of the velocity of the cell phase, unr ; and (f) the bulk pressure, p, throughout the tumour. All plots are at time t = 75. All parameters are as given in Table 3.

Fig. 15. Cross sections of various quantities taken from Fig. 14. (a) Cross section along the line y = 0 and (b) cross-section along the line x = 0 for a non-uniform external applied nutrient.

3416

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

Z

Acknowledgements





$  hk rk



X

 wk þ

X

The authors thank Professor Helen Byrne and Dr. Reuben O’Dea of the School of Mathematical Sciences, University of Nottingham, and Dr. Sarah Waters of Oxford Centre for Industrial and Applied Mathematics, University of Oxford, for helpful discussions on the application of the work presented here. This research was enabled by a studentship from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom awarded through the Life Sciences Interface Doctoral Training Centre at the University of Oxford (Grant No. EP/ E501605/1). Appendix A. Derivation of the weak form of the generalised Stokes problem In this Appendix we provide more details on the derivation of the weak form of the velocity system, Eqs. (27) and (28), which may be written as Eqs. (32) and (33). We begin by deriving some identities that will be required. Suppose r is a symmetric tensor, and w is a vector. On using Eq. (7) for the definition of the divergence of a tensor, Eq. (12) for the definition of the gradient of a vector, and Eq. (20) for the definition of the tensor product of two tensors, we see that

 X @ rij X @ @wi ð$  rÞ  w ¼ wi ¼ ðrij wi Þ  rij @xj @xj @xj i;j i;j  X @ @wi ¼ as rij ¼ rji ðrji wi Þ  rij @xj @xj i;j ¼ $  ðrwÞ  r : $w:

We now use Eqs. (A.1), (A.4) and (A.5) to write this equation as

Z





$  hk rk wk  hk pI þ sk þ Rk þ

X

X

þ p$hk þ

kl

l–k

Z





$  hk rk wk ¼

X

Z



 hk rk wk  ndS ¼

Z

@X

¼

Z



 hk rk n  wk dS

@X



k



hk r n  wk dS þ

@ XDk

Z



 hk rk n  wk dS:

@ XNk

As wk 2 H1Ek0 ; wk ¼ 0 on the Dirichlet boundary oXDk. Applying the boundary condition on the Neumann boundary oXNk, Eq. (14), yields

Z





$  hk rk wk ¼

Z

X

sk  wk dS:

ðA:7Þ

@ XNk

Using Eq. (A.2) we see that

hk ð$wk Þ : sk 

X

ðA:8Þ

ðpðhl $hk  hk $hl Þ þ wkl hl $hk  wlk hk $hl þ bkl hk hl ðul  uk ÞÞ:

l–k

X

bkl hk hl ðul  uk Þ  wk

l–k

 pðhk $  wk þ wk  $hk ÞdV ! Z Z X ¼ sk  wk dS þ hk Rk þ hl wkl $  wk

ðA:2Þ

@ XNk

X

X

l–k

ðwkl hl $hk  wlk hk $hl Þ  wk dV:

l–k

Noting that wkl = wlk, and noting that

hk $  wk þ wk  $hk ¼ $  ðhk wk Þ; allows us to write

Z X

Noting that Eq. (3) implies that, for k = 1, 2, . . ., N,

hl ¼ 1  hk ;

hk ð$wk Þ : sk  m bkl hk hl ðul  uk Þ  wk  p$  ðhk wk ÞdV l–k ! Z Z X ¼ sk  wk dS þ hk Rk þ hl wkl $  wk @ XNk

l–k

þ

we may write Eq. (A.3) as

X

ðwkl hl $hk  wlk hk $hl þ bkl hk hl ðul  uk ÞÞ  wk dV ¼ 0: ðA:6Þ

þ

X

!

Using the fact that rk is symmetric and applying the divergence theorem to the first term in Eq. (A.6) yields

Z

ðA:3Þ X

X

Substituting Eqs. (A.7) and (A.8) into Eq. (A.6) yields, after some manipulation,

To begin our derivation of the weak form of Eq. (6), we first consider the summation term in this equation. On using the definitions of fkl and pkl given by Eqs. (8) and (9), we may write, for k = 1, 2, . . ., N:

f ¼

: $wk

l–k

i;j

X

! !

hl wkl I

I : $wk ¼ $  wk :

ðA:1Þ

I ij rij ¼ traceðrÞ:

X l–k

Another vector identity that we shall require is that, for any tensor r,

I :r¼

kl

f  wk dV ¼ 0:

l–k

X

X

l–k

wkl ðhl $hk  hk $hl Þ  wk dV:

ðA:9Þ

l–k kl

f ¼ pðð1  hk Þ$hk  hk $ð1  hk ÞÞ

l–k

þ

X

ðwkl hl $hk  wlk hk $hl þ bkl hk hl ðul  uk ÞÞ

l–k

¼ p$hk þ

X

ðwkl hl $hk  wlk hk $hl þ bkl hk hl ðul  uk ÞÞ:

ðA:4Þ

l–k

Using Eqs. (10), (11) and (26), we may write the stress tensor for phase k, k = 1, 2, . . ., N, in terms of the overall mixture pressure, p, and viscous component of the stress tensor for phase k, sk, as

rk ¼ pI þ sk þ Rk þ

X

! hl wkl I :

ðA:5Þ

l–k

We now take the scalar product of Eq. (6) with a test function wk 2 H1Ek0 , where the space H1Ek0 is defined in Section 3.3.1, and integrate over the domain X:

This equation is identical to Eqs. (22)–(25), the weak form of Eq. (6). Eq. (A.9) may be summed over all phases k = 1, 2, . . ., N to give the weak formulation given by Eq. (27), as described in Section 3.3.2. We now turn our attention to writing down the weak form of Eq. (21). Multiplying this equation by q 2 L2(X) and integrating over the domain X yields

Z

$

X

N X

! hk uk q dV ¼ 0:

k¼1

Interchanging the order of summation and integration allows us to write N Z X k¼1

q$  ðhk uk Þ dV ¼ 0;

X

from which the weak form given by Eq. (28) follows.

J.M. Osborne, J.P. Whiteley / Comput. Methods Appl. Mech. Engrg. 199 (2010) 3402–3417

References [1] G. Lemon, J.R. King, H.M. Byrne, O.E. Jensen, K.M. Shakesheff, Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory, J. Math. Biol. 52 (2006) 571–594. [2] R.D. O’Dea, S.L. Waters, H.M. Byrne, A two-fluid model for tissue growth within a dynamic flow environment, Eur. J. Appl. Math. 19 (2008) 607–624. [3] R.D. O’Dea, S.L. Waters, H.M. Byrne, A multiphase model for tissue construct growth in a perfusion bioreactor, Math. Med. Biol. 27 (2010) 95–127. [4] B.G. Sengers, C.W.J. Oomens, F.P.T. Baaijens, An integrated finite-element approach to mechanics, transport and biosynthesis in tissue engineering, ASME J. Biomech. Engrg. 126 (2004) 82–91. [5] C.J.W. Breward, H.M. Byrne, C.E. Lewis, The role of cell–cell interactions in a two-phase model for avascular tumour growth, J. Math. Biol. 45 (2002) 125– 152. [6] K.A. Landman, C.P. Please, Tumour dynamics and necrosis: surface tension and stability, IMA J. Math. Appl. Med. 18 (2001) 131–158. [7] C.J.W. Breward, H.M. Byrne, C.E. Lewis, A multiphase model describing vascular tumour growth, Bull. Math. Biol. 65 (2003) 609–640. [8] W. Alt, M. Dembo, Cytoplasm dynamics and cell motion: two-phase flow models, Math. Biosci. 156 (1999) 207–228. [9] J. Jung, A. Hassanein, R.W. Lyczkowski, Hemodynamic computation using multiphase flow dynamics in a right coronary artery, Ann. Biomed. Engrg. 34 (2006) 393–407. [10] V.C. Mow, S.C. Kuei, W.M. Lai, C.G. Armstrong, Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiment, ASME J. Biomech. Engrg. 102 (1980) 73–84.

3417

[11] S.M. Wise, J.S. Lowengrub, H.B. Frieboes, V. Cristini, Three-dimensional multispecies nonlinear growth—I model and numerical method, J. Theor. Biol. 253 (2008) 524–543. [12] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, UK, 1994. [13] J. VonNeumann, R. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys. 21 (1950) 232–237. [14] D.J. Acheson, Elementary Fluid Dynamics, Clarendon Press, Oxford, UK, 1990. [15] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford University Press, Oxford, UK, 2005. [16] N. Bölgen, Y. Yang, P. Korkusuz, E. Güzel, A.J. El Haj, E. Pisßkin, Threedimensional ingrowth of bone cells within biodegradable cryogel scaffolds in bioreactors at different regimes, Tissue Engrg. A 14 (2008) 1743–1750. [17] A.J. El Haj, S.L. Minter, S.C. Rawlinson, R. Suswillo, L.E. Lanyon, Cellular responses to mechanical loading in vitro, J. Bone Miner. Res. 5 (1990) 923–932. [18] A.M. Freyria, Y. Yang, H. Chajra, C.F. Rousseau, M.C. Ronziere, D. Herbage, A.J. El Haj, Optimization of dynamic culture conditions: effects on biosynthetic activities of chondrocytes grown in collagen sponges, Tissue Engrg. 11 (2005) 674–684. [19] J.M. Osborne, R.D. O’Dea, J.P. Whiteley, H.M. Byrne, S.L. Waters, The influence of bioreactor geometry and the mechanical environment on engineered tissues, ASME J. Biomech. Engrg. 132 (2010) 051006. [20] H.M. Byrne, J.R. King, D.L.S. McElwain, L. Preziosi, A two-phase model of solid tumour growth, Appl. Math. Lett. 16 (2003) 567–573. [21] G. Helmlinger, P.A. Netti, H.C. Lichtenbeld, R.J. Melder, R.K. Jain, Solid stress inhibits the growth of multicellular tumour spheroids, Nat. Biotechnol. 5 (1997) 778–783.