The Biot–Stokes coupling using total pressure: Formulation, analysis and application to interfacial flow in the eye

Abstract

We consider a multiphysics model for the flow of Newtonian fluid coupled with Biot consolidation equations through an interface, and incorporating total pressure as an unknown in the poroelastic region. A new mixed-primal finite element scheme is proposed solving for the pairs fluid velocity–pressure and displacement–total poroelastic pressure using Stokes-stable elements, and where the formulation does not require Lagrange multipliers to set up the usual transmission conditions on the interface. The stability and well-posedness of the continuous and semi-discrete problems are analysed in detail. Our numerical study is framed in the context of applicative problems pertaining to heterogeneous geophysical flows and to eye poromechanics. For the latter, we investigate different interfacial flow regimes in Cartesian and axisymmetric coordinates that could eventually help describe early morphologic changes associated with glaucoma development in canine species.

Introduction

Poroelastic structures are found in many applications of industrial and scientific relevance. Examples include the interaction between soft permeable tissue and blood flow, or the study of the spatial growth of biofilm in fluids. When the interaction with a free fluid is considered, the mechanics of the fluid and poroelastic domains are coupled through balance of forces and continuity conditions that adopt diverse forms depending on the expected behaviour in the specific application (see, e.g., [1], [2], [3], [4] and the references therein). The particular problem we consider in this paper as motivation for the design of the finite element formulation is the interfacial flow of aqueous humour between the anterior chamber and the trabecular meshwork (which is a deformable porous structure) in the eye, and how such phenomenon relates to early stages of glaucoma.

Glaucoma encompasses a group of mechanisms that lead to decreased retinal function, impaired visual fields and blindness. The main risk factor for glaucoma in canines is an abnormal increase in the intra-ocular pressure (which under physiologically normal conditions is balanced between aqueous humour production and outflow to the venous drainage system [5]). We are interested in modelling the flow behaviour of aqueous humour within the anterior chamber and its interaction with the poroelastic properties of particular compartments in the drainage outlet located between the base of the iris and the limbus, which, in the dog eye and most other non-primate species consists of an array of thin tissue columns (pectinate ligaments) [6] which mark the boundary of the trabecular meshwork with the anterior chamber. Sketches of the regions of interest are depicted in Fig. 1.1. Our focus is on how the physical changes associated with pectinate ligament dysplasia, a change seen in all dogs with primary angle closure glaucoma, affect aqueous humour flow through this boundary. We stress that the ciliary cleft anatomy of all carnivorous mammals is fairly similar to that of the dog. In fact, this anatomy is preserved across much of herbivorous mammals as well (see, e.g., [7]). Therefore, even though the dog is probably the most studied due to its status as a companion animal, this work likely applies to most carnivorous and herbivorous mammals.

The flow within the anterior chamber will be modelled by Navier–Stokes and Stokes’ law for Newtonian fluids, whereas the filtration of aqueous humour through the deformable trabecular meshwork and towards the angular aqueous plexus will be described by Darcy’s law. Pressure differences are generated by production (from the ciliary muscle) and drainage (to angular aqueous plexus and then linked to the veins at the surface of the sclera through collecting channels) of aqueous humour.

Other effects that could contribute to modification of the flow patterns and that we do not consider here, are thermal properties (buoyancy mechanisms due to temperature gradients from inner to outer cornea) [8], cross-link interaction between fibrils in the cornea [9], pressure changes due to phacodonesis (vibration of the lens while the head or eye itself moves) and Rapid Eye Movement during sleep [10], and nonlinear flow conditions in the filtration region (incorporated in [11] through Darcy–Forchheimer models).

In contrast with [12], [13], [14], [15], here we consider that the coalescing of the pectinate ligaments results in marked changes in porosity properties of the anterior chamber — trabecular meshwork interface, which could eventually lead to progressive collapse of the ciliary cleft. We further postulate that these modifications of the tissue’s microstructure could be induced by forces exerted by the flow that concentrate at the interface between the dysplastic pectinate ligament and the anterior chamber, and which occur over a timescale much larger than that of the ocular pulsating flow. In fact, evidence of the compliance of the trabecular meshwork can be found in, e.g., [16]. One of the earliest modelling works including a coupling between aqueous humour in the anterior chamber with complying structures is presented in [17], where mechanical properties of the bovine iris were employed to set an elastic interface to represent blinking. Other fluid–structure interaction models have been recently developed in [18], suggesting that flow conditions in the trabecular meshwork and the outlets could be largely affected by the changes of permeability in microstructure, and [19], where poroelastic properties of the choroid and viscoelastic response of the vitreous body are used to set up a more complete 3D model of larger scale that discards a dedicated physiological description of the trabecular meshwork and considers instead a windkessel model.

In the general context of single phase fluid / poromechanical coupling, there are already a variety of finite element formulations starting from the work [20], which focuses on the effects of secondary consolidation. More recently, partitioned finite element formulations using domain decomposition and or Nitsche’s approach for single and double poroelastic layers in contact with a single phase fluid can be found in [21], [22], [23], [24]. Monolithic couplings have been analysed in [25] for a mixed Darcy formulation using a Lagrange multiplier to impose flux continuity (see also [26] for the extension to the case of non-Newtonian fluids), as well in [27], [28] for a primal Darcy formulation. Ghost penalty methods have been employed for cut FEM methods valid in the regime of large deformations in [29].

Here, and drawing inspiration from the formulation in [25], [30], [31], [32], we rewrite the poroelasticity equations using three fields (displacement, fluid pressure and total pressure). Compared to previous works [21], [22], [23], [24], [25], [26], [27], [28], which employ the classical displacement formulation, an advantage of the present approach, inherited from [32], is that the formulation is free of poroelastic locking, meaning that it is robust with respect to the Lamé parameters of the poroelastic structure. This is of particular importance when we test variations of the flow response to changes in the material properties of the skeleton and when the solid approaches the incompressibility limit. The present work also stands as an extension of the formulation recently employed in [3] (where only the case of intrinsic incompressible constituent in the poroelastic region were considered) to obtain approximate solutions for heterogeneous poroelasticity coupled with Stokes flow in channels (and using also heterogeneous elastic moduli); while the PDE analysis, numerical aspects, and applicability of the formalism to more realistic scenarios have not yet been addressed. In this work, under adequate assumptions, the analysis of the weak formulation is carried out, using a (time continuous) semi-discrete Galerkin approximation and a weak compactness argument. The well-posedness of the semi-discrete formulation is established using the theory of differential algebraic equations (see, e.g., [33]) and using similar results to those obtained in [25], [34]. A conforming mixed finite element scheme of general order is used. Furthermore, a fully discrete scheme based on backward Euler’s time discretisation is considered, and the unique solvability and convergence for the fully discrete scheme are established.

We have organised the contents of this paper in the following manner. Section 2 outlines the model problem, motivating each term in the balance equations and stating the interfacial and boundary conditions. Section 3 states the weak form of the governing equations in Cartesian and axisymmetric coordinates. Then, in Section 4, we address the construction of the finite element scheme, the well-posedness of the continuous and discrete problems, the stability of the fully discrete system in matrix form. Section 5 states the fully-discrete scheme and presents the error estimates. In Section 6 we collect computational results consisting in verification of spatio-temporal convergence and analysis of different cases on simplified and more physiologically accurate geometries, including also a typical application in reservoir modelling. One of the examples involves large displacements near the interface, in which case a harmonic extension operator is used to deform the fluid domain. We close with a summary, some remarks and a discussion on model generalisations in Section 7.