A Hybrid Reactive Multiphasic Mixture With a Compressible Fluid Solvent

The fields of biomechanics and biophysics increasingly require more sophisticated models that can be used to describe complex physiological processes between molecules, cells, tissues, and organs within living organisms. Mixture theory has been previously used to model mixtures of solid, fluid, and solute constituents [1–4], which had led to significant advances in characterizing the mechanics of biological tissues and cells [5–10]. The biphasic theory, which represents one of the earlier applications of mixture theory to biological tissues, models a binary mixture of an intrinsically incompressible fluid flowing through a porous solid matrix whose skeleton material is intrinsically incompressible [5]. The triphasic theory was extended from biphasic theory [5] by also including the monovalent cation and anion of a single salt [7]. Multiphasic theory was subsequently formulated by generalizing the assumptions of the triphasic theory, incorporating any number of neutral or charged solutes to facilitate the modeling of the mechanoelectrochemical behavior of tissues as well as reactive processes [11–13]. Variations of multiphasic theory have been used in diverse fields such as articular cartilage mechanics [14–17], osmotic loading of cells [18], solute transport induced by dynamic loading of porous tissues [19–21], drug transport in nervous tissue [22,23], intervertebral disk [24] and arterial walls [25], modeling atherosclerosis in arteries [26], and brain swelling and edema [27].

However, there are still some important problems in biomechanics that cannot be easily modeled using existing mixture frameworks. Here, we illustrate two such problem areas, to serve as motivation for extending the existing multiphasic mixture framework: (a) mechanobiological responses of cells subjected to the flow of a biological fluid, such as blood flow over endothelial cells [28], synovial fluid flow over synoviocytes and chondrocytes [29,30], and interstitial fluid flow around osteocytes in bone [31,32], and (b) transport and sedimentation of aerosols and liquid droplets in air flows [33,34], including flow of expired air through deformable porous face masks, or flow of ventilation air through purifying air filters [35–37]. Both of these categories of problems require modeling the flow of a viscous fluid at Reynolds numbers up to $∼104$⁠, where the viscous fluid (blood, synovial fluid, expired air, etc.) may be a mixture of a solvent (water or air) and multiple solutes of various molecular weights, diffusivities, and concentrations (salt ions, cytokines, aerosols, liquid droplets), with Peclet numbers up to $∼108$⁠. The flow may interact with deformable solid or porous structures (biological cells, face mask or air filter), producing significant tractions and deformations, and these interactions may lead to reactive processes in the fluid and the deformable porous or solid domains with which it interacts (stretch-activated, voltage-activated, or osmotically activated cell membrane channels leading to a cascade of intracellular signaling events, liquid droplet evaporation in air, or binding of droplets and aerosols to the mask).

The main impediment against using the existing standard multiphasic theory [11,12], which we previously implemented in the open-source finite element software febio [10,38–40], is its inability to accommodate dynamic flow of a viscous fluid. More recently however, we implemented a computational fluid dynamics solver in febio where the fluid undergoes isothermal processes [41]. To accommodate fluid–structure interactions, we formulated a specialized mixture framework whereby the viscous fluid flows within a porous deformable solid domain with zero fluid–solid frictional interactions; the deformable solid domain exhibits negligible but nonzero stiffness, sufficient to regularize the mesh deformation [42]. Then we extended this theoretical framework to allow dynamic flow of the viscous fluid through a deformable porous solid domain that has non-negligible mass, stiffness, and fluid–solid frictional interactions. As suggested in prior work [43], we called this framework ‘hybrid biphasic theory’ [44], because it treats the fluid as a compressible medium (contrary to the standard biphasic theory [5]), while the skeleton of the porous solid matrix remains intrinsically incompressible. More recently, we reported the results of our finite element implementation of this framework in febio [45]. This hybrid biphasic theory represents the starting point for the formulation of a more general hybrid reactive multiphasic theory as proposed here, where we include solutes and reactive processes to accommodate the range of applications listed in the previous paragraph.

Therefore, the objective of this study is to extend the hybrid biphasic theory to incorporate a charged solid matrix, neutral and charged solutes, and reactions taking place among the various mixture constituents, while also accounting for the dynamics of the solid and viscous solvent constituents. The most significant challenge in this study, in comparison to our prior formulations of the hybrid biphasic theory [44] and standard multiphasic theory [12,39], is the accommodation of osmotic effects when the solvent is compressible. The acceleration of solutes is neglected, since such terms are generally insignificant and their incorporation would require including the velocity of each solute as additional unknowns. Electroneutrality prevails in this formulation and we assume that modeled processes are isothermal, as these assumptions facilitate the formulation of jump conditions for the axioms of mass, momentum and energy balance across interfaces, such as the interface between a fluid solution and a hybrid reactive multiphasic mixture.

The hybrid multiphasic mixture consists of a porous solid matrix and an interstitial fluid that includes a solvent and multiple solutes, undergoing isothermal processes. The skeleton of the porous solid matrix is intrinsically incompressible, but the pores may lose or gain volume. The solvent is modeled as a compressible viscous fluid. Solutes are modeled as intrinsically incompressible inviscid fluid constituents, having negligible volume fraction in the fluid solution. The solvent, typically water, is taken to be electrically neutral and, though it may serve as a catalyst in reactive processes, it does not gain or lose mass from reactions. As a result of these modeling assumptions, a hybrid multiphasic mixture may change in volume when subjected to a hydrostatic pressure, due exclusively to the compressibility of the solvent. As in the case of the related hybrid biphasic mixture [44], waves may propagate at two distinctive speeds through the compressible porous solid matrix and the compressible fluid solution.

where $ρTrf$ is the true solvent density in the solvent's reference configuration, dV^f is the elemental volume of the pore fluid in the current configuration, and $dVrf$ is its referential volume [45].

In these expressions, $ρTσ$ is constant for all σ due to their intrinsic incompressibility.

Remark 1. In this formulation we preclude changes in the solvent mass resulting from reactions, implying that the referential true solvent density, $ρTrf$⁠, always remains constant. In the presence of reactions that alter the concentrations of solid bound molecules σ, such as growth and remodeling processes, $ρTs$ may evolve according to Eq. (2.2) due to reactive processes [39].

Remark 2. By definition, in a constrained mixture of solids all solid constituents σ share the same velocity, $vσ=vs$ [9].

where $DαfDt=∂f∂t+gradf·vα$ is the material time derivative of any arbitrary function $f(x,t)$ in the spatial frame, following the motion of constituent α.

Remark 3. In a solid mixture, individual mixture constituents σ may have a deformation gradient $Fσ$ that differs from $Fs$⁠, as would occur in growth and remodeling processes where solid constituents may get deposited at different times. $Fσ$ may be related to $Fs$ by constitutive assumptions about the growth process, as reviewed in Ref. [46]. The simplest constitutive model is that $Fσ=Fs$ for all σ. Other models can be found in Refs. [47] and [48].

Remark 4. According to the general framework of mixture theory we may define a fluid deformation gradient $Ff=∂χf/∂Xf$⁠, consistent with the definition for the solid. By analogy to the solid, $Ff$ describes the deformation of the porous fluid (whose pores are filled with the solid matrix). Therefore, its determinant would represent the volumetric strain of this porous fluid, which does not have the same meaning as J^f as defined in Eq. (2.1). Thus, the determinant of $Ff$ is not a useful state variable for assessing the intrinsic compressibility of the fluid, so we choose not to use it. Therefore, we allow ourselves to use the variable name J^f in Eq. (2.1) as a measure of fluid compressibility, recognizing that it represents a minor abuse of notation since it may be misconstrued as the determinant of $Ff$⁠.

Thus, $φf$ is entirely determined by $φrs$ (which may evolve due to reactions involving the solid [39]), and the solid volume ratio J^s, which we can calculate from the solid motion.

is the volumetric flux of the fluid relative to the solid [10,12,44]. In this mixture framework we need to satisfy the mixture mass balance (2.9) and all of the solid and solute mass balances.

Remark 5. In the standard reactive multiphasic theory, where the solvent is also modeled as intrinsically incompressible (⁠$Jf=1$⁠), the mixture mass balance reduces from Eq. (2.9) to $div(vs+w)=∑αρ̂α/ρTα$ [39].

where $Vα=Mα/ρTα$ is the molar volume of α. Since solutes are assumed to occupy a negligible volume fraction in the mixture, reactions involving solutes are similarly assumed to cause no significant change in the mixture volume. Therefore, we may choose to evaluate the above summation only over the solid constituents $α=σ$ (since we assume that the solvent is not exchanging mass in reactive processes, $ρ̂f=0$⁠).

where $zα$ is the charge number of constituent α. The first of these equations represents the standard electroneutrality condition, valid in the current configuration, whereas the second enforces electroneutrality even as reactions are taking place. The latter requirement is equivalent to enforcing $∑αzανα=0$ for each reaction.

or equivalently $divIe=0$⁠, where $Ie≡Fc∑ιzιjι$ is the electric current density and F_c is Faraday's constant. In numerical implementations such as the finite element method, this current density constraint must be enforced explicitly in order to produce a stable solution for the electric potential of a mixture.

emphasizing that these supplies represent internal forces that must cancel out [12].

Remark 7. For each solute momentum balance equation the inertia term $ριaι$ will be neglected in the application of this hybrid multiphasic framework to biological tissues.

Remark 8. Consistent with Remark 2, all solid constituents σ in a constrained solid mixture share the same diffusion velocity $uσ=us$⁠. The stress in individual solid constituents is indeterminate due to internal reaction forces that enforce the constraint $vσ=vs$⁠. However, these internal reaction forces cancel out when evaluating the sum of stresses in all the solid constituents [46].

introduced into the entropy inequality (2.28) using the Lagrange multiplier $Fcψ$ [8,12]. As becomes apparent later, ψ represents the electric potential in the mixture.

Remark 9. In the standard multiphasic theory, where all the constituents are assumed to be intrinsically incompressible (implying $Jf=1$ and $ρTf=constant$ in the current treatment), the mixture mass balance obtained from the sum of Eq. (2.4) divided by $ρTα$ over all α reduces to $∑αdiv(φαvα)−ρ̂α/ρTα=0$ when making use of the mixture saturation condition $∑αφα=1$⁠. In its alternative form, $∑α(ραI:Lα+gradρα·uα−ρ̂α)/ρTα=0$⁠, this equation is introduced as a constraint in the entropy inequality, using a Lagrange multiplier p which ends up representing the fluid pressure in the mixture [12].

For our hybrid multiphasic theory we want to accommodate stresses and pressures arising from the solid deformation and fluid volumetric strain, so we include $Fs$ and J^f in our list of state variables. We also account for viscous stresses in the fluid, requiring the addition of $Lf$ to that list. Furthermore, we want to account for evolving solid composition by including $ρrσ$⁠, and evolving solute concentrations, leading us to also add $ρrι=ριJs$ to the list of state variables. Finally, we account for frictional interactions between some of the constituents, requiring us to include the diffusive velocities $uα$⁠. Since we plan to use this framework to examine only isothermal processes, we exclude θ and its spatial gradient from our state variable list, which now reduces to $(Fs,Jf,Lf,ρrσ,ρrι,uα)$⁠. Based on the principle of equipresence, a priori all functions of state in this formulation (i.e., $ρ̂α, σα$⁠, $p̂α, qα, ε̂α$⁠, $ψα$⁠, and $ηα$⁠) depend on this complete list of state variables.

Thus, μ^f appearing in Eq. (2.34) is based on the above relation.

From Eq. (2.39) we conclude that the (inner part of the) referential free energy density $ΨrI$ only depends on the state variables $(Fs,Jf,ρrσ,ρrι)$⁠. Since the entropy is zero, we also conclude that the internal energy is equal to the free energy of the mixture.

Since solutes occupy a negligible fraction of the mixture and are intrinsically incompressible, we now assume that their specific free energies only depend on solute concentrations, $ψβ=ψβ(cι)$⁠, where $β≠s,f$ spans all solutes ι. We also assume that the solid and solvent specific free energies (⁠$ψσ$ and ψ^f) depend on $cι$⁠, as solute concentrations control the osmolarity of the interstitial fluid, whereas charged solutes may also interact with the fixed charge density of the solid constituents (e.g., producing charge shielding).

Remark 10. According to the relations presented above, $ρrι=(Js−φrs)Mιcι$⁠. Here, J^s evolves with $Fs$ and $φrs$ evolves with $ρrσ$⁠, which are state variables independent of $ρrι$⁠. Therefore, we can treat $cι$ as an independent state variable used in lieu of $ρrι$⁠, such that their differential increments satisfy $dρrι=(Js−φrs)Mιdcι$⁠.

Remark 11. The solvent chemical potential μ^f includes a measure of the fluid pressure, $p/ρTf$⁠, consistent with prior standard mixture formulations that employ an intrinsically incompressible solvent [7,8,11,12]. To emphasize the contribution of the pressure p to the chemical potential μ^f, we have previously chosen to refer to this quantity as the mechano-chemical potential [46], by analogy to the concept of electrochemical potential. However, the expression on the right-hand side of Eq. (2.53) has also been described as the specific Gibbs function (or Gibbs potential, or Gibbs energy) in the classical thermodynamic literature, though a more modern terminology calls it the specific free enthalpy. It should be emphasized that all these names refer to the same physical measure given here.

for the solutes, where inertia terms have been neglected.

Constitutive relations for the specific free energies $ψα$ appearing in Eq. (2.49), the mass density supplies $ρ̂α$ in reactive processes, the viscous stress $τf$⁠, and the dissipative momentum supplies $p̂dα$⁠, may be application-dependent. However, they must satisfy the residual dissipation statement in Eq. (2.43). The first two terms of this dissipation statement, $τf:Lf−∑αp̂dα·uα$⁠, represent the heat dissipation resulting from frictional interactions within the solvent (⁠$τf$⁠) and friction between constituents (⁠$p̂dα$⁠). The remaining terms in the residual dissipation statement determine if the reactions associated with the mass density supplies $ρ̂α$ can take place spontaneously (i.e., without the addition of heat, $∑αρ̂α(μα+12uα·uα)+ρTfμf∑σρ̂σ/ρTσ≤0$⁠), or if they require additional heat, e.g., as supplied by the frictional dissipative mechanisms in Eq. (2.43). Other forms of heat supply may be modeled in a more general framework that allows spatio-temporal variation of temperature. Regardless of the presence of reactive processes, we may formulate constitutive relations for $τf$ and $p̂dα$ that enforce $τf:Lf≥0$ and $∑αp̂dα·uα≤0$⁠, which are sufficient (but not necessary) conditions to satisfy the residual dissipation statement (2.43).

This relation is equivalent to Onsager's reciprocity relation in the classical field of irreversible thermodynamics [7]. To satisfy the residual dissipation statement uniquely, it suffices to let $fαβ$ be symmetric and positive semidefinite.

Remark 12. For many applications of mixture theory to porous-permeable biological soft tissues, it can be shown that $τf:Lf≪−∑αp̂dα·uα$ in the residual dissipation statement (2.43). Therefore, in those cases it is common to prescribe $τf≈0$⁠.

In these expressions, R is the universal gas constant and θ is the absolute temperature. All remaining $fαβ$'s for frictional interactions between solutes are taken to be zero. While each of the functions of state k, $d0ι$⁠, and $dι$ may depend on any of the state variables $(Fs,Jf,Lf,ρrσ,cι,uα)$⁠, in biological tissues k varies primarily with the solid matrix strain (thus, a strain measure associated with $Fs$⁠) [56–59], $dι$ varies with the strain and solute concentrations, and $d0ι$ varies with solute concentration [12,53]. More generally, k, $d0ι$ and $dι$ may also be modeled such that they evolve with mixture composition.

These last two relations represent the standard physical chemistry expressions for the chemical potential of solutes and solvent in an ideal solution.