The effect of magnetic field on flow induced-deformation in absorbing porous tissues

1.   Introduction

When fluid flows through a rigid porous material, the dynamics is adequately described either by Darcy's law ^[1] for low speed flow or by nonlinear Forcheimer's law ^[2] for high speed flow. On the other hand, for flow through a deformable porous material a more sophisticated model is required to study the interaction between elastic solid and fluid. In fact, the properties such as permeability and porosity of these materials alter due to the forces associated with the flow, which in turn, affect the passage of fluid through these materials. The pioneering work on deformation of a porous material coupled with fluid flow dates back to Terzaghi ^[3] in 1925. Later on, it was extended by Biot to study the soil consolidation ^[4] and soil culmination ^[5] and these works have been used extensively in soil mechanics.

Following the pioneering work of Fick ^[6] on mixture theory, Truesdell ^[7] put forth the modern theory by proposing balance equations appropriate to mixtures irrespective of their constitution. This theory was then comprehensively reviewed and applied to study chemically inert mixture of two ideal gases by Atkin and Craine ^[8] along with applications ^[9]. Later on, Rajagopal and Tao ^[10] authored an excellent book on mechanics of mixtures in which conservation laws were derived and several example problems were discussed from mixture theory viewpoint. In addition to these, a review of mixture theory for deformable porous media along with applications in different scientific fields was presented very recently by Siddique et al. ^[11].

The application of mixture theory to soft biological tissues essentially started with the work of Kenyon ^[12,13], who discussed radial flux of fluid through a porous cylinder replicating a model of flow through arterial tissue. Following this, Jayaraman ^[14] studied the problem of flow through an artery wall with constant permeability, and Jain and Jayaraman ^[15] also investigated the same problem but considered two layers in the artery wall each with different permeabilities. Similarly, Klanchar and Tarbell ^[16] studied the water flow through arterial tissue by considering a linear form of the strain-dependent permeability. Apart from arterial tissue, the theory was extended further by Mow et al. ^[17,18,19] by considering the articular cartilage as a deformable porous material saturated with synovial fluid. The main application, however, was the lubrication properties of synovial joints such as knee. Due to compression of articular cartilage, the synovial fluid moves through the pores of the tissue to form a surface lubricating layer. Very recently, Ahmed et al. ^[20] developed a mathematical model of non-Newtonian flow-induced deformation in absorbing porous tissues using mixture theory. Other biological tissues that have been modeled in this way include the skin ^[21], the lung ^[22] and the cornea ^[23], etc.

In clinical medicine the prediction of absorption rate of a substance from an injection site is a problem of fundamental and great importance. Therefore, understanding the mechanical process involved in this mechanism may help resolve many unclear and unresolved issues. Although, we consider here a model with certain simplifying assumptions, nevertheless, this model can be used as a precursor to models of a number of other clinically important problems. The problem of fluid absorption from an injection site into the brain tissue was investigated by Nicholson ^[24] by considering two possibilities for fluid flow through the tissue. Firstly, that the injected fluid saturates a region of the tissue and then perfuses to the remaining tissue and secondly that a small cavity of injected fluid forms initially followed by the fluid flow through rest of the tissue. We consider here the second possibility to model the growth of a spherical cavity with time under a time dependent pressure applied within the cavity.

In magnetohydrodynamic (MHD) processes, the motions of an electrically conducting fluid induce and maintain a finite magnetic field. In particular, these MHD processes have become relevant to many applications in medical science such as cell isolation, drug targeting, tissue activation and magnetization, magnetic hyperthermia and magnetic resonance imaging. In the past few years, a magnetohydrodynamic theory for deformable porous media has gained a considerable attention of many researchers. Eldabe et al. ^[25] investigated the magnetohydrodynamic flow of a bi-viscosity fluid through a porous medium in a layer of deformable material. In particular, the authors used a mixture theory approach to develop a model for flow through a symmetrical channel. Following this idea, several authors ^[26,27] studied the problem of infiltration of MHD liquid into a deformable porous material. Sreenadh et al. ^[28] investigated MHD Couette flow of a Jeffrey fluid over a deformable porous layer in a channel bounded below by a finite deformable porous layer and by a moving rigid plate using a two component mixture theory. The applications of magnetic field to biological systems has also attracted the attention of many authors in the last few decades. To this end, application of magnetic fields in artificial human joints were developed through the models incorporating bone-in growth in porous implants ^[29,30]. Moreover, the externally applied magnetic field stimulates the functions of biological tissues along with their regeneration ^[31]. Apart from these studies, the effect of magnetic field on synovial joints was investigated by Tandon et al. ^[32] using lubrication theory. Specifically, the authors considered the knee joint as two approaching porous cartilageous surfaces and the synovial fluid in the cartilage was represented by an Oldroyd model. It was reported that a suitably designed magnetic field improves the performance of synovial joints along with a better articulation in a diseased state.

The aim of this study is to develop a mathematical model of flow-induced deformation from pressurized cavities in absorbing porous tissues subject to a uniform applied magnetic field. In particular, a model with a spherical cavity embedded in a porous medium of infinite extent is used to find fluid pressure and solid displacement of the tissue as a function of radial distance and time. The governing set of equations are nondimensionalized using suitable dimensionless quantities which are then solved numerically to assess the influence of magnetic parameter. The arrangement of this article is as follows: in section 2, we present mathematical formulation of the problem. Section 3 deals with the solution methodology for the governing set of equations followed by results and discussion in section 4. Finally, concluding remarks are presented in section 5.

2.   Mathematical formulation

In order to develop the model, we use a biphasic mixture theory to describe the fluid flow and solid deformation in a biological tissue under the influence of a uniform magnetic field as shown in Figure 1. The main underlying idea of mixture theory is that each constituent of the mixture is continuous and occupies every point in space at each instant of time ^[8]. The orientation of the applied magnetic field is taken in such a way that it only produces the radial flow in the tissue from the injection site.

We assume that injected fluid in the porous tissue is conducting and viscous while the organic solid matrix is isotropic, homogeneous and linearly elastic. The gravitational and osmotic forces are neglected. It is also assumed that shear stresses are negligible on account of one dimensional radial flow and the constituents of the mixture are intrinsically incompressible. Under these assumptions, equations of motion for solid and fluid phase in the radial direction are written as (see Appendix for details).

where $\phi^s$ and $\phi^\ell$ represent solid and fluid volume fractions and $v^s$ and $v^\ell$ are velocities of solid and fluid phase, respectively, $\rho^\ell_T$ is the intrinsic density of fluid phase, $\beta$ is the proportionality constant which depends upon the concentration of capillaries and lymphatics in the tissue and permeability of their walls, $p$ is the fluid pressure, $K$ the drag coefficient of relative motion, $\sigma_{0}$ the electric conductivity of the fluid, $B_0$ the uniform magnetic flux. It is important to note that the last term on right hand side of equation (2.3) is the contribution of MHD and left hand side of equation (2.4) is the divergence of solid stress in the radial direction where $\sigma_{rr}$ and $\sigma_{\theta \theta} = \sigma_{\phi \phi}$ are defined as components of solid stress ^[33]. Note that equations (2.1) and (2.2) represent the mass balance relations for solid and fluid phase, respectively, whereas equations (2.3) and (2.4) are derived from solid and fluid momentum balances ^[25,26]. The term appearing on right hand side of equation (2.2) is due to loss of fluid mass at a rate proportional to fluid pressure while it passes through capillaries and lymphatics. Adding equations (2.1) and (2.2), we obtain

which on denoting the macroscopic medium velocity in the radial direction, $v_r = \phi^{s} v^s+ \phi^{\ell} v^{\ell}$, and using the relation, $\phi^s+\phi^\ell = 1$, reduces to

Using equation (2.4) into (2.3) and simplifying, yields

Denoting $v^s = \frac{\partial u}{\partial t}$, where $u$ is a component of solid displacement and using an expression for $v^\ell$ from macroscopic medium velocity $v_r$ into the relation (2.7) and keeping in view equation (2.4), we have

where

is defined to be the permeability of the medium ^[33]. Equation (2.8) can be explained physically by taking into account the Darcy's law and considering the solid stress as being governed by the standard equilibrium equation of the theory of linear elasticity. As suggested by the relation ($\ref{ten eq}$), the permeability $\kappa$ of the porous medium decreases as a consequence of solid compression. Taking into account equation (2.4) and the relation $\frac{\partial \sigma_{rr}}{\partial r} + 2 \frac{\sigma_{rr}-\sigma_{\theta \theta}}{r} = H_{a} \frac{\partial \phi}{\partial r}$ where $H_{a} = \lambda+2 \mu$ is the agammaegate modulus and $(\lambda, \mu)$ are Lamé constants (see Appendix for details), the governing equation (2.8) can be written as

which on equating the first two expressions and then integrating gives the relationship between $p$ and $\phi$ as

where both $p$ and $\phi$ tend to zero as $r\rightarrow\infty$. From equations (2.10) and (2.11), a relation for $v_r$ in terms of fluid pressure $p$ and solid displacement $u$ may be obtained as

Combining equations (2.6) and (2.12), we obtain

As the porous tissue is compressed, consequently its permeability would reduce due to a decrease in porosity. The relationship between porosity and displacement of the medium is given by ^[33]

where $\bf G$ is the left Cauchy-Green deformation tensor and $\phi^s_{0}$ is the initial solid volume fraction. For infinitesimal solid deformations this can be approximated as ^[33,34]

where $\phi = \nabla \cdot \bf{u}$ is the one-dimensional dilatation and $\alpha_{0} = \frac{\phi^s_{0}}{\phi^{\ell}_{0}}$. The permeability as a function of porosity can then be expressed as $\kappa = \kappa(\phi)$ which is related to fluid pressure via equation (2.11). Recognizing the first term on left hand side of equation (2.13) as the time derivative of $\phi = \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 u)$ and keeping in view the relation (2.11), we obtain

This equation is expressed in terms of fluid pressure $p$ except the third term on left hand side which still involves a component of solid displacement $u$ in the parenthesis. In order to eliminate this component, we in view of relations (2.11) and (2.15) assume a linear permeability of the form

where $k_0$ and $n$ are material constants. Note that this simple form of permeability, which is a general approximation to any $\kappa(\phi)$ for small $\phi$, is valid for infinitesimal solid deformations and various authors ^[16,35] have considered a similar form to study the water flow through arterial tissue and radial flow through deformable porous shells. Using the relations (2.11), (2.15) and (2.17) into the governing equation (2.16), we obtain

Assuming initial fluid volume fraction $\phi^{\ell}_{0}$ to be constant and $n = \alpha_{0}$, this equation can be reduced to a more convenient form as

This equation may be non-dimensionalized using the following dimensionless quantities

where $t_0, r_0, p_0$ and $u_0$ are typical time, radius, pressure and deformation scales, respectively. After introducing these choices, equation (2.19) on dropping the bars takes the following form

where various dimensionless parameters are defined as

Note that values of various parameters appearing in this equation may differ considerably among biological tissues ^[33]. The equation (2.21) is required governing equation in terms of non-dimensional fluid pressure $p(r, t)$ under the action of a uniform applied magnetic field and is related to porosity $\phi$ via equation (2.11). This parabolic PDE is nonlinear which makes it difficult to solve analytically. We thus use a numerical method for its solution which will be described in the next section. It is worth mentioning that setting the magnetic parameter $M = 0$ in equation (2.21) and assuming a linear permeability of the form (2.17), we recover the Newtonian fluid case ^[33]. A natural time scale $t_0$ for the current problem may be obtained by setting $\frac{\alpha}{1+M} = 1$ as

which reduces the parameter $\omega$ in equation (2.22) to

The non-dimensional boundary conditions for fluid pressure $p(r, t)$ and solid displacement $u(r, t)$ are written as

and

where $f(r)$ and $g(t)$ are some specified functions of their arguments and $L = \frac{\lambda}{\lambda+2 \mu}$. Note that the given displacement boundary condition for infinitesimal deformation was derived from a general boundary condition between a fluid and a porous material ^[36]. Since the governing equation (2.21) is solved for fluid pressure, therefore an equation is required that relates solid displacement to fluid pressure. This is accomplished by combining the relations $\phi = \frac{1}{r^2} \frac{\partial}{\partial r} (r^2 u)$ and $p(r, t) = H_{a} \phi(r, t)$ with equation (2.20) to give

where $\gamma = \frac{p_0 r_0}{H_a u_0}$ is a dimensionless parameter. Integrating this equation and applying the boundary conditions (2.26) for solid displacement, yields

which gives displacement of solid as a function of non-dimensional radial distance $r$ and time $t$. Note that the effect of magnetic parameter on steady-state solutions is absent as suggested by equation (2.21). In order to carry out numerical computations for the nonlinear problem, we assume that $a(t) = 1, g(t) = 1, f(r) = 0, \alpha = 1, \gamma = 1, L = 0.5$ and $\delta = 1$, which are consistent with the values considered by Barry and Aldis ^[33].

3.   Solution methodology

In this section, we briefly outline the procedure for the solution of governing set of equations. To begin with, we first give an exact solution of equation (2.21) in terms of MHD fluid pressure $p(r, t)$ subject to boundary conditions (2.25) for the constant permeability case (i.e. when $n = 0$) and compare the result with method of lines (MOL) as a validation of numerical scheme to be employed for the solution of nonlinear problem. Following the method adopted by Barry and Aldis ^[33] for linearized problem, the solution for fluid pressure for the choice of applied pressure, $g(t) = e^{\frac{-\omega}{1+M} t} t^{\frac{m}{2}}$, is given as

where $m$ is a non-negative integer, $\Gamma(\cdot)$ is the gamma function and ${\bf{i}}$erfc represents integration of the error function ^[37]. As stated earlier, for the nonlinear problem we use method of lines whose main idea is to discretize the space variable and its derivatives and leaving the time variable continuous ^[38]. This space discretization results into a system of coupled ODEs which may then be solved using Matlab's well established and efficient solvers as an initial value problem. For completeness, we now show a brief calculation of the MOL by discretizing the governing equation (2.21) using central finite difference formulas for the first and the second space derivatives as

where

Here $p_0$ represents the left boundary condition and $p_{N+1}$ denotes the right boundary condition of fluid pressure. From the initial condition for fluid pressure, we have

We thus have an initial value problem consisting of N ODEs in equation (3.2) and corresponding initial conditions outlined in equation (3.4) which may now be solved using the efficient solver such as ode23s. It is important to note that in order to carry out numerical simulations, we truncate the spatial domain for fluid pressure at $r = 13$ due to satisfaction of the far field boundary condition and the temporal domain at $t = 7$ because after this value there is no significant change in the fluid pressure and hence the long time behavior of the solution may be obtained. Moreover, solid displacement $u(r, t)$ in equation (2.28) is computed using trapezoidal rule by solving a definite integral involving fluid pressure. It is also important to note that in the implementation of numerical method for the nonlinear problem the left boundary condition was applied at non-dimensional radius $r = 1$ corresponding to initial radius.

In Figure 2 on the left, we present a comparison between the exact (i.e. equation (3.1)) and numerical (i.e. equation (3.2)) solution for the non-dimensional fluid pressure in the presence of magnetic field for the constant permeability case. An excellent agreement between the two solutions can be observed validating the proposed numerical scheme. On the other hand, the graph on the right presents the mesh convergence study for the fluid pressure profile along with the grid independence test. It can be seen that the solution converges rapidly as the number of spatial nodes $N$ increases and after $N = 80$ the solution is not affected by changing the mesh size.

4.   Results and discussion

This section contains the outcome of our numerical simulations for fluid pressure and solid displacement for various values of magnetic parameter. Methods described in the previous section are used and the results are illustrated graphically.

In Figure 3, non-dimensional fluid pressure $p(r, t)$ is plotted against time $t$ for the choice of four different applied pressure profiles $g(t)$ in the presence of magnetic effects. The solid curve on the graph indicates that fluid pressure in the tissue drops off exponentially with time whereas the dashed curve suggests that there is an initial increase in the pressure to a local maximum followed by an exponential decay. Moreover, the dashed-dotted curve on this graph representing a periodic form of applied pressure shows that fluid pressure in the porous tissue rises initially to a maximum value followed by a periodic decrease while the dotted curve representing a quadratic type of applied pressure indicates that pressure of fluid in the tissue is increasing from the center of the cavity. It is to be noted that a similar form of periodic and quadratic pressure was also considered in ^[39].

Figure 4 describes the influence of magnetic parameter $M$ on fluid pressure $p(r, t)$ at $t = 0.5, 7.0$ for a given distance from center of the cavity. This graph shows that fluid pressure in the tissue drops off more rapidly as strength of the magnetic field increases. The Lorentz force associated with applied magnetic field boosts the fluid flow in the tissue which consequently reduces fluid pressure in the porous material. Thus, application of a suitably designed magnetic field may possibly assist to control the fluid flow in deformable porous tissues for practical purposes. On the other hand, for fixed $M$, fluid pressure in the tissue rises with time and there is not much difference in pressure distribution after $t\geq7$ explaining the long time behavior of the solution. Since fluid pressure and porosity of the tissue are related directly via equation (2.11), so this plot also illustrates the porosity as a function of radial distance and time.

The effect of magnetic parameter $M$ on solid displacement $u(r, t)$ at $t = 0.5, 7.0$ is presented in Figure 5. The relation for solid displacement (2.28) suggests that higher fluid pressure in the tissue should induce greater solid deformation, and vice versa. Accordingly, due to a reduction in fluid pressure with magnetic parameter, the tissue solid displacement also decays with magnetic parameter as shown in Figure 5. This means that a properly applied high magnetic field can prevent large solid deformation of the material. Additionally, this graph also illustrates that deformation of the solid increases with time, and that after $t\geq7$, there is not significant change in the displacement distribution of the tissue. Moreover, in the absence of magnetic effects (i.e. $M = 0$) and for large $n$, it is expected that solid displacement $u$ would exhibit greater inflection due to nonlinearities in the system, particularly, for large time.

Finally, non-dimensional pressure-dependent permeability, $\kappa(p) = 1+n \delta p$, is plotted against radial distance $r$ in Figure 6 for different values of magnetic parameter $M$ when $n = 0.25$ and $\delta = 1$. The tissue permeability reduces as strength of magnetic field increases validating the direct relation with fluid pressure. This fact indicates that a high magnetic field would allow less fluid to seep through the porous tissue.

5.   Concluding remarks

In this study we analyzed a one-dimensional model of flow-induced deformation from pressurized cavities in absorbing porous biological tissues under an applied magnetic field. The tissue was assumed to be isotropic, homogeneous and linearly elastic. We used a biphasic mixture theory to develop the model with an assumption of infinitesimal solid deformation. A linear permeability relation was considered which allowed the governing equation to be written explicitly in terms of fluid pressure. A method of lines approach was adopted to solve the nonlinear parabolic PDE in terms of fluid pressure which was then used to find tissue displacement by employing the trapezoidal rule.

We noticed a reasonable reduction in fluid pressure, solid displacement and permeability of the tissue due to the presence of magnetic effects. Thus, a properly designed magnetic field may help achieve the required fluid flow and solid deformation in the tissue for certain physiological applications. The solid deformation enhances fluid absorption in the tissue and alters the porosity and permeability of the material. Although, we assumed that the absorption rate and tissue pressure are related linearly, however, more sophisticated forms of the absorption process may be considered to enhance the understanding of these complex systems.

This work is relevant to the interpretation of experimental studies on neuropharmacology and in situ electrochemistry especially with MHD effects. Other applications include the activation and magnetization of a variety of soft tissues for clinical purposes.

This research serves as a starting point in MHD modeling of flow-induced deformation in absorbing porous tissues and we hope that this simple mathematical model can be further extended in different directions. These directions include consideration of nonlinear stress-strain relation, different permeability relations, variable magnetic field and more realistic geometries. Moreover, incorporation of material isotropy and chemical effects are the topics of further research.

Acknowledgments

The author J. I. Siddique would like to acknowledge the support from Simons Foundation Grant No. 281839.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

We consider a deformable porous tissue as a continuous binary mixture of intrinsic incompressible solid and fluid, where each point in the mixture is occupied by both fluid and isotropic solid. The balance of mass for solid and fluid phase can be written as ^[33]

where ${\bf v}^s$ and ${\bf v}^\ell$ are velocities and ${\rho^s}$ and ${\rho^\ell}$ are densities of the solid and fluid phase, respectively, $p$ is the fluid pressure and $\beta$ is a proportionality constant. The conservation of linear momentum for the $\eta$ phase in the absence of body forces except the Lorentz force is written as ^[25]

where $\eta = s, \ell$ represents either the solid or the fluid phase, ${{\bf T}^{\eta}}$ is the stress tensor for the $\eta$ phase, ${\mathit{\boldsymbol{\pi}}^{\eta}}$ is the friction force between the mixture constituents which satisfies the relation $\mathit{\boldsymbol{\pi}}^s+\mathit{\boldsymbol{\pi}}^{\ell} = {\bf 0}$. Note that ${\bf J }\times {\bf B}$ represents the Lorentz force due to applied magnetic field in which ${\bf{J}}$ designates the current density and ${\bf{B}}$ denotes the magnetic flux density. The stress tensor ${{\bf T}^{\eta}}$ and the drag force ${\mathit{\boldsymbol{\pi}}^{\eta}}$ in equation (A.3) is defined as ^[33]

where ${\bf I}$ is the identity tensor, $K$ the drag coefficient of relative motion, and $\phi^{\eta}$ and $\mathit{\boldsymbol{\sigma}}^{\eta}$ represent the volume fraction and the stress for the $\eta$ phase, respectively. The Maxwell's equations of electromagnetism along with the Ohm's law are written as ^[25,26]

where $\mu_{c}$ is the permeability of free space, $\bf{E}$ the electric field and $\sigma_{0}$ the electric conductivity of the fluid. The term ${\bf{J \times B}}$ in the momentum equation (A.3) can be written in view of Ohm's law as

where the total magnetic field ${\bf{B}}$ may be decomposed as, ${\bf{B}} = {\bf{B_0}+{\bf{b}}}$, in which ${\bf{B_0}}$ is the imposed magnetic field and b is the induced magnetic field which may be ignored on account of low magnetic field Reynolds number approximation. Thus, equation (A.7) when induced magnetic and electric fields are negligible, takes the form

Application of the vector identity $({\bf{X}} \times {\bf{Y}}) \times {\bf{Z}} = {\bf{Y}}({\bf{X}} \cdot {\bf{Z}})-{\bf{X}}({\bf{Y}}\cdot{\bf{Z}})$ reduces equation (A.8) to the following form

which on assuming ${\bf{v}}^{\eta} \cdot {\bf{B_0}} = 0$ leads to

where $B_0$ is the strength of the applied magnetic field ${\bf{B_0}}$. Thus, the momentum equation (A.3) after neglecting the inertial terms and taking into account the relation (A.10) takes the following form

In order to derive the equations (2.3) and (2.4), we proceed as follows. Using the relation (A.4) into equation (A.11) and adding both phase equations, yields

where the relation $\phi^s+\phi^\ell = 1$ has been used. Note that while deriving equation (A.12), the viscous stress ${\mathit{\boldsymbol{\sigma}}^{\ell}}$ is assumed negligible on account of one-dimensional radial flow and ${\mathit{\boldsymbol{\sigma}}^{s}} = \mathit{\boldsymbol{\sigma}}$ is considered for the rest of the derivation ^[20]. Equation (A.12) under the assumptions $\mathit{\boldsymbol{\pi}}^s+\mathit{\boldsymbol{\pi}}^{\ell} = {\bf 0}$ and ${\bf{v}}^{s}+{\bf{v}}^{\ell} = {\bf 0}$ enables us to write

which yields equation (2.4) in the scalar form upon using the definition of the divergence in the radial direction. Substituting the relations (A.4) and (A.5) into equation (A.11) and performing some mathematical manipulation, we have respectively for the solid and fluid phase

Eliminating the pressure $p$ from these equations and using the formula $\phi^s+\phi^\ell = 1$, we obtain

which on using again the assumption ${\bf{v}}^{s}+{\bf{v}}^{\ell} = {\bf 0}$, yields

This equation in combination with the relation (A.13) gives equation (2.3) in component form.

We now turn our attention to the motion of the solid phase. Under the assumptions of infinitesimal deformations and one-dimensional radial flow, the components of stress for the solid phase are defined as ^[33]

where $\lambda$ and $\mu$ are Lamé stress constants and due to spherical symmetry of the problem all other stress components are assumed to be zero. The divergence of stress in the radial direction is given by

Substituting equations (A.18) and (A.19) into equation (A.20) and simplifying, we obtain

where

is the local change in porosity and $H_{a} = \lambda+2 \mu$ is the agammaegate modulus.