Discontinuous Galerkin Methodsfor Mass Transfer through Semi-Permeable Membranes

# Discontinuous Galerkin Methods for Mass Transfer through Semi-Permeable Membranes

Andrea Cangiani Department of Mathematics, University of Leicester, University Road, Leicester LE1 7RH, United Kingdom, e-mail: andrea.cangiani@le.ac.uk    Emmanuil H. Georgoulis Department of Mathematics, University of Leicester, University Road, Leicester LE1 7RH, United Kingdom, e-mail: Emmanuil.Georgoulis@le.ac.uk    Max Jensen Department of Mathematical Sciences, University of Durham, Science Laboratories, South Road, Durham DH1 3LE, United Kingdom, e-mail: m.p.j.jensen@durham.ac.uk
###### Abstract

A discontinuous Galerkin (dG) method for the numerical solution of initial/boundary value multi-compartment partial differential equation (PDE) models, interconnected with interface conditions, is presented and analysed. The study of interface problems is motivated by models of mass transfer of solutes through semi-permeable membranes. More specifically, a model problem consisting of a system of semilinear parabolic advection-diffusion-reaction partial differential equations in each compartment, equipped with respective initial and boundary conditions, is considered. Nonlinear interface conditions modelling selective permeability, congestion and partial reflection are applied to the compartment interfaces. An interior penalty dG method is presented for this problem and it is analysed in the space-discrete setting. The a priori analysis shows that the method yields optimal a priori bounds, provided the exact solution is sufficiently smooth. Numerical experiments indicate agreement with the theoretical bounds and highlight the stability of the numerical method in the advection-dominated regime.

I

nterface modelling, mass transfer, discontinuous Galerkin methods, semilinear parabolic problems, nonlinear interface conditions

## 1 Introduction

Models of mass transfer of substances (solutes) through semi-permeable membranes appear in various contexts, such as biomedical and chemical engineering applications [26]. Examples include the modelling of electrokinetic flows, solute dynamics across arterial walls, and cellular signal transduction (see, e.g., [10, 53, 15] and the references therein).

This work is concerned with the development and analysis of numerical methods for a class of continuum models for mass transfer based on initial/boundary value multi-compartment partial differential equation (PDE) problems, closed by nonlinear interface conditions. The interface conditions considered are the Kedem-Katchalsky (KK) equations, which represent an established model for the mass transfer mechanisms [37, 36]. More specifically, we consider a generic model problem consisting of a system of semilinear advection-diffusion-reaction parabolic PDE problems in multi-compartment configurations, coupled with nonlinear interface KK-type conditions. The focus is to address some challenges in the numerical solution of these models, such as the treatment of nonlinearities due to both the interface modelling and the nonlinear reactions, the discontinuity of the state variables across the interface, as well as the development of stable numerical methods in the advection-dominated regime.

Numerical methods for mass transfer problems based on conforming finite elements have been developed for the solution of solute dynamics across arterial walls; see [53, 46, 45] and the references therein for more details. Some existence results for the purely diffusing interface problem without forcing, coupled with KK-type interface conditions, along with some numerical experiments are given in [14]. Further, numerical approaches to the treatment of interface conditions for PDE problems, resulting to globally continuous solutions can be found, e.g., in [6, 3, 16, 44, 41].

Discontinuous Galerkin (dG) methods (see, e.g., [31, 47, 18] and the references therein) and mortar methods (see, e.g., [9] for a survey) have also been proposed for the treatment of coupled systems via interface conditions in various contexts [50, 27, 28, 22, 23, 33, 30, 29]. Also, the advantages of dG methods for interfacing different numerical methods (numerical interfaces) have been identified [43, 17], as well as their use on transmission-type/high-contrast problems, yielding continuous solutions across the transmission interface, has been investigated [21, 12, 24, 13].

Here, we consider a dG method of interior penalty type for the solution of the semilinear parabolic advection-diffusion-reaction PDE system coupled with nonlinear interface conditions of KK-type across the subdomains. The use of dG is motivated partly by the observation that the interface conditions, yielding discontinuous solutions across the interface, can be imposed by modifying the interior penalty dG numerical fluxes. Another important factor for employing a dG method is the desired stability property of the numerical method in the advection-dominated regime. A priori bounds for the proposed spatially discrete dG method in both the - and -type norms are presented for a range of reaction fields, under the simplifying assumption that the finite element mesh is aligned with the subdomain interfaces.

A priori error bounds for interior penalty dG methods for parabolic problems have been considered in various settings (see, e.g., [47] for an exposition and the more recent [19]). DG methods for semilinear parabolic spatially self-adjoint problems with locally Lipschitz continuous nonlinearity have been analysed in [40]. In the present analysis, advection terms are included and systems of equations are considered. In the presence of advection, the analysis of the symmetric interior penalty dG method in [40] would require the assumption of quasi-uniformity of the mesh. To avoid this assumption, a different continuation argument is employed in the derivation of the a priori bounds presented here, at the expense of a stricter growth condition on the nonlinearity of the forcing term. This continuation argument is inspired by the derivation of a posteriori bounds for semilinear parabolic phase-field models [38, 8]. The nonlinear interface terms are tackled using a non-standard elliptic projection which is inspired by a classical construction of Douglas and Dupont [20] for the treatment of nonlinear boundary conditions.

The paper is organised as follows. In Section 2, the PDE model along with a short derivation of the nonlinear interface conditions is presented. Section 3 is devoted to the description of the dG method proposed for the advection-diffusion part of the spatial operator incorporating the nonlinear interface conditions. Section 4 contains error estimates for the new elliptic projection, which is, in turn, utilised in the subsequent a priori error analysis presented in Section 5. Section 6 contains numerical experiments highlighting the stability and the optimal rate of convergence of the proposed method in practice. Finally, some conclusions are offered in Section 7.

## 2 Interface modelling and governing PDEs

We shall consider systems of parabolic semilinear PDEs describing the flux of solutes around and through a semi-permeable membrane. The membrane is modelled as an internal boundary equipped with nonlinear interface conditions which are described in the following section.

### 2.1 Interface modelling

We outline the Kedem-Katchalsky (KK) equations modelling solutes flow across semi-permeable membranes. The KK equations have been introduced in [37]; we refer to [51] for earlier works and to [26] for a thorough exposition.

It is assumed that the membrane separates two compartments and filled with a free fluid which is called the solvent and that the membrane characteristics are uniform in space and time. The KK equations specify the dependence of the solutes and solvent fluxes across the membrane in terms of two driving forces, namely the hydrostatic and osmotic pressure jumps. In the case of a single solute, the solvent and solute fluxes from to normal to the membrane walls are given by, respectively,

 Jv = LP(δp−σdδπ), (1) Js = ωδπ+Jv(1−σf)¯u, (2)

in terms of the hydrostatic and osmotic pressure jumps and , with being the ideal gas constant, denoting temperature, and the solute concentration jump across the membrane. Here, represents the average concentration of the solute across the membrane. The above constitutive laws are characterised by the phenomenological coefficients of filtration , reflection and , and permeation . These coefficients may depend on the concentration while they are assumed to be constant with respect to both the position along the membrane and time.

Equation (1), which is known as Starling’s law of filtration, shows that the solvent flow is affected by the osmotic flow of the solute. This observation introduces a nonlinearity in the transport term of the solute flux in equation (2). Indeed, substituting into (2) we get the final model for the solute flux:

 Js=ωRTδu+LP(δp−σdRTδu)(1−σf)¯u=p(u1,u2)δu−r¯u(\boldmathb⋅\boldmathn)|Ω2.

In this last expression, we have collected the diffusive part of the flux by introducing the nonlinear permeability function , and written the advective transport in terms of the friction coefficient and the normal component of the transport field yielded by the hydrostatic pressure, denoted by .

In presence of solutes, it is usually assumed that there is no direct coupling between the solutes’ fluxes [26]. Under this simplifying assumption the KK equations describing the flux of the solvent and solutes read

 Jv = LP(δp−n∑j=1σd,jδπj), (3) Js,i = ωδπi+Jv(1−σf)¯ui,i=1,…,n, (4)

with , , and representing, respectively, the flux, osmotic pressure, and average concentration of the -th solute. Proceeding as in the case of one solute we get the following expression for the solute fluxes:

 Js,i = pi(u1i,u2i)δui−n∑j=1~pi,j(\boldmathu1,\boldmathu2)δuj−ri¯ui(\boldmathb⋅\boldmathn)|Ω2 (5) = \boldmathpi(\boldmathu1,% \boldmathu2)⋅δ\boldmathu−ri¯ui(\boldmathb⋅\boldmathn)|Ω2,i=1,…,n,

with , , , where denotes the concentration of the -th solute in the -th compartment, the corresponding permeability, and the cross-coefficients expressing the crowding effect inside the membrane. For the second equality, we have collected in the vector permeability function all diffusion terms.

It remains to fix a model for the average concentrations inside the membrane. In the case of relatively thick membranes [53, 45], it is appropriate to consider a mechanical approach describing solutes flow within the membrane through an advection-diffusion model. This approach yields a weighted arithmetic average of the concentration at the membrane’s faces. Thus, denoting by and , , such concentrations, we get

 ¯ui=δwui:=v1iu1i+v2iu2i, (6)

for some given weights with . These can be expressed in terms of the internal Pèclet number of the membrane advection-diffusion model and are so that the upwind value dominates [53], cf. the conditions given below in (13).

In what follows, the fluxes given by (5) together with the model (6) for the average concentration inside the membrane are used to close the PDE problem with appropriate interface conditions.

### 2.2 Notation

We denote by , , the standard Lebesgue spaces, , , with corresponding norms ; if , we shall write instead . Also the norm of will be denoted by and if by for brevity; by we write the standard -inner product on ; when the arguments are vectors of -functions, the -inner product is modified in the standard fashion. We denote by the standard Hilbertian Sobolev space of index of real-valued functions defined on ; in particular signifies the space of functions in whose traces onto the boundary vanish. For , we denote the standard Bochner spaces , with being a Banach space with norm . Finally, we denote by the space of continuous functions with norm .

Let be a bounded open domain with Lipschitz boundary in , and let be the boundary of . The domain is subdivided into two subdomains and , such that , where , see Figure 1. For , we assume that has Lipschitz boundary and that has positive -dimensional (Hausdorff) measure. We define , .

We shall employ the following notational convention: vectors are indicated with lower case bold symbols, diagonal matrices with upper case (non-bold) symbols, and tensors with upper case bold symbols.

The gradient of a vector function in is a mapping gained from componentwise application of the gradient operation: . Similarly the divergence of the tensor-valued function is where the are rows of .

### 2.3 Model Problem

For a time interval , , and for

 \boldmathu∈L2(0,T;H1),\boldmath% ut∈L2(0,T;H−1),

with , we consider the system of semilinear parabolic equations

 \boldmathut−∇⋅(A∇\boldmathu−U\boldmathB)+\boldmathF(\boldmathu)=\boldmath0in (0,T]×(Ω1∪Ω2), (7)

with denoting the diagonal matrix . Here, is an tensor field with rows , , and diagonal, with , where , . We assume that there exists a constant of uniform parabolicity such that for all and . For simplicity, we also require that the matrix is positive semi-definite. Finally, is a vector field satisfying the growth condition

 |\boldmathF(\boldmathw)−\boldmathF(\boldmathv)|≤C(1+|\boldmathw|+|% \boldmathv|)γ|\boldmathw−\boldmathv|, (8)

for and constant, where denotes the Euclidean distance on . The admissible values of the constant will be discussed in detail at various instances in the text.

We impose the initial condition

 \boldmathu(0,x)=\boldmathu0(x) on {0}×Ω, (9)

for . On , we impose mixed Dirichlet and Neumann boundary conditions as follows. Let be an index running over the set . For each the boundary is split into , with being of positive -dimensional (Hausdorff) measure. Further, we subdivide , where and are the inflow and outflow parts of the boundary for the -th equation. Finally, we assign

 ui =g\rm Dion Γ\rm Di, (10) ai∇ui⋅\boldmathn =g\rm Nion Γ\rm Ni∩∂Ω+i, (ai∇ui−BTiui)⋅%\boldmath$n$ =g\rm Nion Γ\rm Ni∩∂Ω−i,

for Dirichlet and Neumann data , , respectively; here and in what follows denotes the unit outward normal vector to . We denote by the characteristic function of . Upon defining , and , we can write (10) collectively as

 \boldmathu=\boldmathg\rm D% on Γ\rm Dand(A∇\boldmathu−X−U\boldmathB)\boldmathn=\boldmathg\rm Non Γ\rm N% , (11)

where and .

The model problem is completed imposing, across , the fluxes described in Section 2.1. In view of (5) and (6), we define the friction coefficients and weights and the permeabilities , , as functions of the traces of from both sides of the interface. The interface conditions read

 (ai∇ui−uiBTi)⋅\boldmathn|Ω1 =\boldmathpi(\boldmathu1,\boldmathu2)⋅(\boldmathu2−% \boldmathu1)−ri(υ1iu1i+υ2iu2i)(Bi\boldmathn)|Ω1,on ΓI, (ai∇ui−uiBTi)⋅\boldmathn|Ω2 =\boldmathpi(\boldmathu1,\boldmathu2)⋅(\boldmathu1−% \boldmathu2)−ri(υ1iu1i+υ2iu2i)(Bi\boldmathn)|Ω2,on ΓI,

where , . Introducing, , , and , the interface conditions can be written in vector notation as

 (A∇\boldmathu−U\boldmathB)\boldmathn|Ω1 =\boldmathP(\boldmathu)(% \boldmathu2−\boldmathu1)−R(Υ1U1+Υ2U2)(\boldmathB\boldmathn)|Ω1,on ΓI, (12) (A∇\boldmathu−U\boldmathB)\boldmathn|Ω2 =\boldmathP(\boldmathu)(% \boldmathu1−\boldmathu2)−R(Υ1U1+Υ2U2)(\boldmathB\boldmathn)|Ω2,on ΓI.

We make the following assumptions on the weights and permeabilities in accordance with Section 2.1. For every , the weights satisfy, for any ,

 υ1i(\boldmathx)+υ2i(\boldmathx)=1,{υ1i(% \boldmathx)≥υ2i(\boldmathx)% if (Bi\boldmathn|∂Ω1)(\boldmathx)≥0,υ1i(\boldmathx)<υ2i(\boldmathx)otherwise. (13)

We let denote the function describing the diffusive flux across the interface , that is

 ~\boldmathp(\boldmathx1,% \boldmathx2):=\boldmathP(\boldmathx1,\boldmathx2)(\boldmathx1−% \boldmathx2)∀\boldmathx1,% \boldmathx2∈Rn, (14)

and assume that and that its Jacobian is bounded.

Throughout this work, we assume that the problem given by (7), (9), (11), and (12) has a unique solution that remains bounded up to, and including, the final time .

## 3 Space discretisation by the discontinuous Galerkin method

### 3.1 Finite element spaces

Let be a shape-regular and locally quasi-uniform subdivision of into disjoint open elements , such that , the skeleton. Further we decompose into three disjoint subsets , where . We assume that the subdivision is constructed via mappings , where are smooth maps with non-singular Jacobian, and is the reference -dimensional simplex or the reference -dimensional (hyper)cube. It is assumed that the union of the closures of the elements forms a covering of the closure of ; i.e., .

For we denote by the set of polynomials of total degree at most if is the reference simplex, and the set of all tensor-product polynomials on of degree in each variable, if is the reference hypercube. Let be given for each . We consider the -discontinuous finite element space

 Vh:={v∈L2(Ω):v|κ∘Fκ∈Pmκ(^κ),κ∈T}, (15)

and set .

Next, we introduce relevant trace operators. Let , be two elements sharing an edge . Denote the outward normal unit vectors on of and by and , respectively. For functions and that may be discontinuous across , we define the following quantities: for , and , on the restriction to , we set , , and , , where denotes the standard tensor product operator, with . If , these definitions are modified as follows: and .

Further, we introduce the mesh quantities , with , , if , and the averaged values , , if . Finally, we define and .

We shall assume the existence of a constant independent of such that, on any face that is not contained in , given the two elements , sharing that face, the diffusion matrix satisfies

 C−1A≤∥∥A∥∞,κ∥A−1∥∥∞,κ′≤CA. (16)

This assumption can be removed using the ideas from [24], but we refrain from doing so here for simplicity of the presentation.

The next result is a modification of the classical trace estimate for functions in ; see [11] for similar results.

{lemma}

Assume that the mesh is both shape-regular and locally quasi-uniform. Then for , the following trace estimate holds:

 2∑j=1∥v|Ωj∥2ΓI≤c1ϵ(∑κ∈T∥∇v∥2κ+∥%\boldmath$h$−1/2[v]∥2Γ\rm int)+c2ϵ−1∥v∥2, (17)

for any and for some constants and , depending only on the shape-regularity of the mesh and on the domain . {proof} We use the decomposition of into a conforming part and a non-conforming part . This decomposition is described in [34, 35] for functions in ; the extension to follows by taking . Using Theorem 2.1(iii) from [35], there exists a , , such that

 ∑κ∈T∥Dα(v−vic)∥2κ∩Ωi≤C∥\boldmathh1/2−|α|[v]∥2Γ%int∩Ωi, (18)

where is the differentiation operator for a multi-index with . Hence, (18) implies

 ∑κ∈T∥Dα(v−vc)∥2κ≤C∥\boldmathh1/2−|α|[v]∥2Γ\rm int, (19)

for , .

The triangle inequality implies

 2∑j=1∥v|Ωj∥2ΓI≤22∑j=1(∥vc|Ωj∥2ΓI+∥(v−vc)|Ωj∥2ΓI). (20)

To bound the first term on the right-hand side of (20), we note that , giving

 2∑j=1∥vc|Ωj∥2ΓI ≤C2∑j=1(∥vc∥2Ωj+∥vc∥Ωj∥∇vc∥Ωj)1/2≤C(∥vc∥2+∥vc∥∥∇vc∥)1/2 (21) ≤C(ϵ−1∥vc∥2+ϵ∥∇vc∥2)1/2,

for any sufficiently small. To further bound the right-hand side of (21), we use the triangle inequality for each term, viz.,

 ∥vc∥≤∥v−vc∥+∥v∥and∥∇vc∥2≤2∑κ∈T(∥∇(v−vc)∥2κ+∥∇v∥2κ), (22)

in conjunction with (19), to arrive at

 2∑j=1∥vc|Ωj∥2ΓI ≤C(ϵ−1∥\boldmathh1/2[v]∥2Γ\rm int+ϵ−1∥v∥2+ϵ∥% \boldmathh−1/2[v]∥2Γ\rm int+ϵ∑κ∈T∥∇v∥2κ)1/2 (23) ≤C(ϵ−1∥v∥2+ϵ∥\boldmathh−1/2[v]∥2Γ\rm int+ϵ∑κ∈T∥∇v∥2κ)1/2,

using the assumption that .

To bound the second term on the right-hand side of (20), we use the trace estimate on each element in conjunction with (19), viz.,

 2∑j=1∥(v−vc)|Ωj∥2ΓI ≤C∑κ∈T:¯κ∩ΓI≠∅(∥\boldmathh−1/2(v−vc)∥2κ+∥\boldmathh1/2∇(v−vc)∥2κ)≤C∥[v]∥2Γ\rm int% . (24)

Noting the assumption , the use of (23) and (24) in (20) concludes the proof.

### 3.2 Space discretization

The discretization of the space variables will be based on a discontinuous Galerkin method of interior penalty type for the diffusion part and of upwind type for the advection. Special care has to be given to the incorporation of the interface conditions.

More specifically, we shall introduce a interior penalty discontinuous Galerkin (IPDG, for short) discretisation of the advection-diffusion operator

 −∇⋅(A∇\boldmathw−W\boldmathB), (25)

where , for .

To this end, we define to be the IPDG bilinear form

 ∑κ∈T∫κ(A∇\boldmathwh−Wh\boldmathB):∇\boldmathvh+∫ΓI({Wh\boldmathB}+BI[[\boldmathwh]]):[[\boldmathvh]] (26) −∫Γ\rm int({A∇\boldmathwh−Wh\boldmathB}:[[\boldmathvh]]+{A∇\boldmathvh}:[[\boldmathwh]]−(Σ+B)[[\boldmathwh]]:[[\boldmathvh]]) −∫Γ\rm D((A∇\boldmathwh−X+Wh\boldmathB):(\boldmathvh⊗\boldmathn)+(A∇\boldmathvh):(\boldmathwh⊗\boldmathn)−Σ\boldmathwh⋅\boldmathvh) +∫Γ\rm N(X+Wh\boldmath% B):(\boldmathvh⊗\boldmathn),

and define

 N(\boldmathwh,\boldmathvh):=∫ΓI(\boldmathP(\boldmathwh)[[\boldmathwh]]−(I−R)({Wh% \boldmathB}+BI[[\boldmathwh]])):[[\boldmathvh]]. (27)

Here, denotes the discontinuity-penalization parameter matrix with constant. Furthermore, and is diagonal with non negative entries.

###### Remark \thetheorem

We comment on the interface terms. The diffusion term, appearing in , is simply given by . This resembles the typical jump stabilisation term with the permeability coefficient replacing the discontinuity-penalization parameter, rendering its implementation within a dG computer code straightforward.

To ensure the coercivity of , the advective interface term has been split as , resulting into contributions in both and . Indeed, the advective interface contribution in can be recast using the weighted mean , so that

 {Wh\boldmathB}υ:[[\boldmathvh]]=({Wh\boldmathB}+BI[[\boldmathwh]]):[[\boldmathvh]], (28)

thereby resembling the typical dG upwinding for linear advection problems and, hence, ensuring the coercivity of .

###### Remark \thetheorem

In this setting, can have non-trivial intersection with both and , thereby extending the discontinuous Galerkin method proposed in [32].

## 4 Elliptic projection error

The a priori error analysis is based on a (non-standard) elliptic projection inspired by a construction of Douglas and Dupont for the treatment of nonlinear boundary conditions in the context of conforming finite element methods [20].

{definition}

For each we define the elliptic projection to be the solution of the problem: find , such that

 (29)

for some fixed , and denoting the exact solution.

The constant in the definition above will be chosen large enough to ensure the uniqueness of the projection , cf. Lemma 4 below.

Next, denoting by , , we define the dG-norm on

 |∥\boldmathw|∥:=( ∑κ∈T(∥√A∇% \boldmathw∥2κ+12∥√\diag(∇⋅%\boldmath$B$)\boldmathw∥2κ)+∥√Σ[[\boldmathw]]∥2Γ\rm D∪Γ\rm int (30) +∥√B[[\boldmathw]]∥2Γ∖ΓI+∥√BI[[\boldmathw]]∥2ΓI)1/2,

where