## Abstract

Mixture theory is a general framework that has been used to model mixtures of solid, fluid, and solute constituents, leading to significant advances in modeling the mechanics of biological tissues and cells. Though versatile and applicable to a wide range of problems in biomechanics and biophysics, standard multiphasic mixture frameworks incorporate neither dynamics of viscous fluids nor fluid compressibility, both of which facilitate the finite element implementation of computational fluid dynamics solvers. This study formulates governing equations for reactive multiphasic mixtures where the interstitial fluid has a solvent which is viscous and compressible. This hybrid reactive multiphasic framework uses state variables that include the deformation gradient of the porous solid matrix, the volumetric strain and rate of deformation of the solvent, the solute concentrations, and the relative velocities between the various constituents. Unlike standard formulations which employ a Lagrange multiplier to model fluid pressure, this framework requires the formulation of a function of state for the pressure, which depends on solvent volumetric strain and solute concentrations. Under isothermal conditions the formulation shows that the solvent volumetric strain remains continuous across interfaces between hybrid multiphasic domains. Apart from the Lagrange multiplier-state function distinction for the fluid pressure, and the ability to accommodate viscous fluid dynamics, this hybrid multiphasic framework remains fully consistent with standard multiphasic formulations previously employed in biomechanics. With these additional features, the hybrid multiphasic mixture theory makes it possible to address a wider range of problems that are important in biomechanics and mechanobiology.

## 1 Introduction

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 $\u223c104$, 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 $\u223c108$. 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.

## 2 Formulation of the Hybrid 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.

### 2.1 Mixture Composition.

*α*, with $\alpha =s$ for the solid, $\alpha =f$ for the fluid solvent, and $\alpha =\iota $ for the solutes. The definitions of the apparent density $\rho \alpha $, the volume fraction $\phi \alpha $, and true density $\rho T\alpha $ for each constituent

*α*are given in prior treatments [12,44] and not repeated here. Solutes occupy a negligible volume fraction, $\phi \iota \u226a1$ and their true densities $\rho T\iota $ are constant, due to their intrinsic incompressibility. The mixture saturation condition (the absence of voids in the mixture) implies that $\phi s+\phi f\u22481$. Since the solvent is compressible, we identify a kinematic state variable that takes into account its volumetric strain. As done in the hybrid biphasic mixture [44,45], we use the solvent volume ratio

*J*relating the solvent current and reference states

^{f}where $\rho 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].

*s*may itself consist of a constrained reactive mixture of intrinsically incompressible solid constituents

*σ*[46], such that its apparent and true mass densities, and its volume fraction are given by

In these expressions, $\rho T\sigma $ 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, $\rho Trf$, always remains constant. In the presence of reactions that alter the concentrations of solid bound molecules *σ*, such as growth and remodeling processes, $\rho 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\sigma =vs$ [9].

### 2.2 Motion and Mass Balance.

where $D\alpha fDt=\u2202f\u2202t+gradf\xb7v\alpha $ 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\sigma $ that differs from $Fs$, as would occur in growth and remodeling processes where solid constituents may get deposited at different times. $F\sigma $ may be related to $Fs$ by constitutive assumptions about the growth process, as reviewed in Ref. [46]. The simplest constitutive model is that $F\sigma =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=\u2202\chi f/\u2202Xf$, 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*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$.

^{f}*α*due to reactions with all other constituents. The temporal evolution of the referential solid mass densities $\rho r\sigma =Js\rho \sigma $ and associated solid volume fractions $\phi r\sigma =Js\phi \sigma $ can been reduced from Eqs. (2.4) and (2.3) as described previously [39]. As in the hybrid biphasic mixture [44], the mixture saturation condition produces

Thus, $\phi f$ is entirely determined by $\phi 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.

*ι*, a constant value. Here, we have also provided a similar expression for the molar concentration supply $c\u0302\iota $. In addition, it is also convenient to define the molar solute flux relative to the solid

*α*, and making use of the mixture saturation condition $\u2211\alpha \phi \alpha =1$, produces the mixture mass balance equation

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)=\u2211\alpha \rho \u0302\alpha /\rho T\alpha $ [39].

*Remark 6*. In a chemical reaction, the mass density supplies $\rho \u0302\alpha $ may be related to the molar production rate $\zeta \u0302$ via $\rho \u0302\alpha =\phi fM\alpha \nu \alpha \zeta \u0302$, where $\nu \alpha $ is the net stoichiometric coefficient (product minus reactant) of

*α*in the reaction [39]. Therefore we may rewrite

where $V\alpha =M\alpha /\rho T\alpha $ 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 $\alpha =\sigma $ (since we assume that the solvent is not exchanging mass in reactive processes, $\rho \u0302f=0$).

#### 2.2.1 Electroneutrality Condition.

where $z\alpha $ 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 $\u2211\alpha z\alpha \nu \alpha =0$ for each reaction.

*c*to the referential fixed charge density $crF$ via

^{F}or equivalently $divIe=0$, where $Ie\u2261Fc\u2211\iota z\iota j\iota $ 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.

### 2.3 Momentum Balance.

*α*, $b\alpha $ is the body force per mass acting on constituent

*α*, $p\u0302\alpha $ is the momentum supply per volume to

*α*due to momentum exchanges with other mixture constituents, and $a\alpha $ is the acceleration of

*α*, given by the material time derivative of $v\alpha $ in the spatial frame, following the motion of

*α*. As shown previously [12,50], applying the axiom of mixtures to obtain the mixture momentum balance produces a mixture stress

*inner part of the mixture*stress, and

*diffusive velocity*of constituent

*α*relative to the mixture. The definition of the mixture velocity

**v**can be found in earlier treatments [1,3]; it represents the mass-weighted barycentric velocity of the mixture constituents. The momentum supplies must satisfy the constraint

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

*Remark 7*. For each solute momentum balance equation the inertia term $\rho \iota a\iota $ 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\sigma =us$. The stress in individual solid constituents is indeterminate due to internal reaction forces that enforce the constraint $v\sigma =vs$. However, these internal reaction forces cancel out when evaluating the sum of stresses in all the solid constituents [46].

### 2.4 Energy Balance.

*α*is repeated here [12,49]

*α*, $L\alpha =gradv\alpha $ is the velocity gradient, $q\alpha $ is the heat flux through

*α*, $r\alpha $ is the specific heat supply to

*α*from sources that are not explicitly modeled, and $\epsilon \u0302\alpha $ is the energy density supply to

*α*due to exchanges with other mixture constituents. As with the axiom of momentum balance, the sum of the individual constituent energy balances (2.22) is equated to the energy balance of a pure substance according to the axiom of mixtures, producing a constraint for the energy density supplies in the form

### 2.5 Entropy Inequality.

*θ*[5,7,44,51]. The mixture entropy inequality is given in the form of the Clausius-Duhem inequality as

*α*. The free energy density of

*α*is then defined by $\Psi \alpha =\rho \alpha \psi \alpha $ where the inner part of the mixture free energy density is given by $\Psi I=\u2211\alpha \Psi \alpha $. The referential free energy density of the mixture is

*α*as $H\alpha =\rho \alpha \eta \alpha $ and the mixture entropy density as $H=\u2211\alpha H\alpha $, where the referential mixture entropy density is

introduced into the entropy inequality (2.28) using the Lagrange multiplier $Fc\psi $ [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 $\rho Tf=constant$ in the current treatment), the mixture mass balance obtained from the sum of Eq. (2.4) divided by $\rho T\alpha $ over all *α* reduces to $\u2211\alpha div(\phi \alpha v\alpha )\u2212\rho \u0302\alpha /\rho T\alpha =0$ when making use of the mixture saturation condition $\u2211\alpha \phi \alpha =1$. In its alternative form, $\u2211\alpha (\rho \alpha I:L\alpha +grad\rho \alpha \xb7u\alpha \u2212\rho \u0302\alpha )/\rho T\alpha =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].

#### 2.5.1 Thermodynamic Restrictions.

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 $\rho r\sigma $, and evolving solute concentrations, leading us to also add $\rho r\iota =\rho \iota 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\alpha $. 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,\rho r\sigma ,\rho r\iota ,u\alpha )$. Based on the principle of equipresence, a priori all functions of state in this formulation (i.e., $\rho \u0302\alpha ,\u2009\sigma \alpha $, $p\u0302\alpha ,\u2009q\alpha ,\u2009\epsilon \u0302\alpha $, $\psi \alpha $, and $\eta \alpha $) depend on this complete list of state variables.

*J*in Eq. (2.1), along with Eq. (2.3), and the mass balance for the solid constituents, the fluid mass balance Eq. (2.4) may be rewritten as

^{f}*α*, while its electrochemical potential is given by

*α*is electrically neutral, $z\alpha =0$ and $\mu \u0303\alpha =\mu \alpha $. Using the definition of Eq. (2.35) for $\alpha =f$, along with Eq. (2.1) and $\rho rf=Js\rho f$, we also note that the chemical potential of the fluid solvent is related to $\u2202\Psi rI/\u2202Jf$ via

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 $\Psi rI$ only depends on the state variables $(Fs,Jf,\rho r\sigma ,\rho r\iota )$. Since the entropy is zero, we also conclude that the internal energy is equal to the free energy of the mixture.

#### 2.5.2 Constitutive Assumptions for 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, $\psi \beta =\psi \beta (c\iota )$, where $\beta \u2260s,f$ spans all solutes *ι*. We also assume that the solid and solvent specific free energies ($\psi \sigma $ and *ψ ^{f}*) depend on $c\iota $, 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, $\rho r\iota =(Js\u2212\phi rs)M\iota c\iota $. Here, *J ^{s}* evolves with $Fs$ and $\phi rs$ evolves with $\rho r\sigma $, which are state variables independent of $\rho r\iota $. Therefore, we can treat $c\iota $ as an independent state variable used in lieu of $\rho r\iota $, such that their differential increments satisfy $d\rho r\iota =(Js\u2212\phi rs)M\iota dc\iota $.

*ψ*depends on

^{f}*J*. Similarly, since the porous solid is deformable, we assume that $\psi \sigma $ depends on $Fs$. Finally, we adopt the simplifying assumption that $\psi \sigma $ does not depend on any of the (potentially evolving) solid concentrations $\rho r\sigma $, implying that solid composition influences the referential free energy density of each solid constituent

^{f}*σ*only based on the proportionality relation $\Psi r\sigma =\rho r\sigma \psi \sigma $. With these constitutive assumptions, it follows from Eq. (2.26) that

*Remark 11*. The solvent chemical potential *μ ^{f}* includes a measure of the fluid pressure, $p/\rho 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

*μ*, we have previously chosen to refer to this quantity as the

^{f}*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.

*β*takes the form

**0**) produces

for the solutes, where inertia terms have been neglected.

### 2.6 Constitutive Relations.

Constitutive relations for the specific free energies $\psi \alpha $ appearing in Eq. (2.49), the mass density supplies $\rho \u0302\alpha $ in reactive processes, the viscous stress $\tau f$, and the dissipative momentum supplies $p\u0302d\alpha $, may be application-dependent. However, they must satisfy the residual dissipation statement in Eq. (2.43). The first two terms of this dissipation statement, $\tau f:Lf\u2212\u2211\alpha p\u0302d\alpha \xb7u\alpha $, represent the heat dissipation resulting from frictional interactions within the solvent ($\tau f$) and friction between constituents ($p\u0302d\alpha $). The remaining terms in the residual dissipation statement determine if the reactions associated with the mass density supplies $\rho \u0302\alpha $ can take place spontaneously (i.e., without the addition of heat, $\u2211\alpha \rho \u0302\alpha (\mu \alpha +12u\alpha \xb7u\alpha )+\rho Tf\mu f\u2211\sigma \rho \u0302\sigma /\rho T\sigma \u22640$), 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 $\tau f$ and $p\u0302d\alpha $ that enforce $\tau f:Lf\u22650$ and $\u2211\alpha p\u0302d\alpha \xb7u\alpha \u22640$, which are sufficient (but not necessary) conditions to satisfy the residual dissipation statement (2.43).

#### 2.6.1 Dissipative Momentum Supplies.

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\alpha \beta $ be symmetric and positive semidefinite.

*Remark 12*. For many applications of mixture theory to porous-permeable biological soft tissues, it can be shown that $\tau f:Lf\u226a\u2212\u2211\alpha p\u0302d\alpha \xb7u\alpha $ in the residual dissipation statement (2.43). Therefore, in those cases it is common to prescribe $\tau f\u22480$.

**k**of the solvent through the porous solid matrix, the solute diffusivity $d0\iota $ in free solution, and the (hindered) solute diffusivity tensor $d\iota $ in the porous-permeable mixture according to

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

#### 2.6.2 Chemical Potential.

*Ideal solutions*: In physical chemistry, for ideal solutions (in the absence of a solid matrix), the free energy of solute

*ι*is given by

*ι*in the solution. Similarly, the free energy of the solvent is given by

*x*is the mole fraction of the solvent. In a solution, mole fractions satisfy $\u2211\alpha \u2260sx\alpha =1$. Thus, the mole fraction of the solvent may also be evaluated as $xf=1\u2212\u2211\iota x\iota $, where the summation is taken over all solutes. For dilute solutions, it follows from a series expansion of the natural logarithm in Eq. (2.69) that

^{f}*ψ*when the solution contains no solutes. Thus, $\psi \iota $ depends only on its own $x\iota $, whereas

^{f}*ψ*depends on all $x\iota $'s. Note that

^{f}*x*is not an independent state variable in this context. For a dilute solution, we can express the mole fraction of solute

^{f}*ι*in terms of its molar concentration $c\iota $ as

*c*

_{0}is an arbitrary constant defining the solute standard state (typically, $c0=1\u2009M$) and $\mu 0\iota $ is a constant that represents the value of $\psi \iota $ when $c\iota =c0$. Similarly, the solvent specific free energy in Eq. (2.70) may be rewritten as

*μ*in Eq. (2.53), we get

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

##### Nonideal Solutions and Mixtures.

*activity*of solute

*ι*and $\mu 0\iota $ is the chemical potential in the standard state when $a\iota =1$. For nonideal solutions, we can express $a\iota $ as

*activity coefficient*$\gamma \iota (Fs,Jf,\rho r\sigma ,c\beta )$ of solute

*ι*is a nondimensional property that describes the deviation of the solute chemical potential from the ideal dilute solution [60] (here,

*β*spans all

*ι*). In Eq. (2.78), $\kappa \iota $ is the

*solubility*of solute

*ι*in the mixture, which may represent the fraction of the interstitial pore volume accessible to the solute [53], such that $0\u2264\kappa \iota \u22641$. For ideal mixtures, we let $\gamma \iota =1$ and for perfect solubility we let $\kappa \iota =1$. Since $\gamma \iota $ and $\kappa \iota $ always appear as a ratio in these equations, we may define the effective solubility [61] as

*ι*may be written as

*effective concentration*is defined as

*osmotic coefficient*of the solvent, a nondimensional property that describes the deviation of the osmotic pressure from the ideal behavior known as van't Hoff's law [64]. From Eq. (2.83), we can define the

*effective pressure*, which represents the “mechanical” or “hydraulic” contribution to

*p*as

whereas $R\theta \Phi \u2211\iota c\iota $ represents the “chemical” or “osmotic” contribution to the fluid pressure.

#### 2.6.3 Implication for Momentum Balances.

### 2.7 Jump Conditions Across Reactive Interfaces.

*α*takes the form

where $u\Gamma \alpha =v\alpha \u2212v\Gamma $, $v\Gamma $ is the velocity of interface Γ, **n** is a unit normal on Γ, and $\rho \xaf\alpha $ is the mass flux of *α* produced by reactions on Γ. In this notation, for any arbitrary vector **v**, we have $[[v]]\xb7n=v+\xb7n++v\u2212\xb7n\u2212=(v+\u2212v\u2212)\xb7n$, such that + and – designate the value of any quantity on the two sides of Γ, $n+=n$ is the outward normal to the + side of Γ, and $n\u2212=\u2212n$. In multiphasic mixtures, we may define the interface Γ on a porous solid boundary, where $v\Gamma =vs$ and $u\Gamma s=0$.

where $c\xaf\iota =\rho \xaf\iota /M\iota $ is the reactive molar flux of *ι* on Γ (per mixture area). This shows that the solute molar flux normal to Γ exhibits a jump whenever reactions take place on Γ which release ($c\xaf\iota >0$) or consume ($c\xaf\iota <0$) that solute. For example, this type of reactive process could be used to model cell membrane channels or transporters that carry solute *ι* across the cell membrane.

**n**in Eq. (2.90) with the unit vector

**s**tangent to Γ [44], producing

From this statement, we can conclude that the traction vector from the inner part of the mixture stress is continuous across Γ.

*ι*is not involved in a reaction on Γ (thus, $\rho \xaf\iota =0$ in Eq. (2.88)). In that case, Eq. (2.94) may be simplified to

*sufficient*(but not necessary) jump conditions are satisfied:

where the solute electrochemical potential is also continuous across Γ, when solute *ι* is not involved in a reaction on Γ. The first two conditions apply when the solvent is present on both sides of Γ, whereas the last condition is applicable when solute *ι* is similarly present on both sides. If solute *ι* is in fact involved in a reaction on Γ, we use Eq. (2.89) to prescribe $c\xaf\iota $ as a flux (Neumann) boundary condition on *ι*, instead of satisfying (2.98) as a Dirichlet boundary condition.

*ψ*that does not depend on solute concentrations. Substituting this expression into the general relation of Eq. (2.50) produces

^{f}This result demonstrates that *μ ^{f}* is exclusively a function of

*J*when the mixture is an ideal solution. It follows from Eq. (2.76) that the effective pressure $p\u0303=p\u2212R\theta \u2211\iota c\iota $ of an ideal solution is also only a function of

^{f}*J*, even though the actual fluid pressure

^{f}*p*depends on

*J*and the solute concentrations, as per (2.52).

^{f}*J*. Based on the above remark regarding ideal solutions, and the fact that $p\u0303$ is the nonosmotic part of the pressure

^{f}*p*as given in Eq. (2.84), it follows that $p\u0303(Jf)$ for nonideal mixtures must also depend only on the fluid volume ratio, implying that $\mu f=\mu f(Jf)$ for nonideal mixtures as well. For example, in this isothermal framework, we may choose a familiar constitutive relation valid for arbitrary (but strictly positive) values of

*J*

^{f}*J*, the continuity of the solvent chemical potential

^{f}*μ*across Γ, expressed in Eq. (2.97), implies continuity of

^{f}*J*, on the understanding that the same solvent species is present on both sides of Γ, thus

^{f}As reported in our prior hybrid biphasic theory [44], the jump condition in Eq. (2.107) suggests that the viscous stress (and thus, the viscosity) of a fluid flowing in a porous solid matrix may scale with the porosity of that medium, such that $\tau f=\phi f\tau $ where $\tau $ would be the viscous stress in the free fluid (absent a porous matrix, when $\phi f=1$). Thus, we could use standard constitutive relations for the viscous stress $\tau $ of Newtonian or non-Newtonian fluids and adapt those models to a multiphasic mixture where $\tau f$ is evaluated from $\phi f\tau $.

## 3 Canonical Problems

The governing equations for reactive hybrid multiphasic media may be applied to a very broad range of problems, most of which require numerical solution methods. In the two canonical problems presented here, we illustrate phenomena fundamentally related to the two representative problem areas described in the introduction. First, we consider the sedimentation of a solute under the action of gravity, as a problem representative of saliva droplets sedimenting in air. Second, we consider the steady Couette flow of a viscous electrolyte in a channel over a charged deformable porous matrix, representing the flow of a biological fluid over the glycocalix of endothelial cells.

*p*, where we insert the constitutive model of Eq. (2.101) for $p\u0303$, thus reducing the solvent chemical potential to the form given in Eq. (2.103). We employ (2.1) to solve for

*ρ*and Eq. (2.6) to solve for $\rho \iota $. When present, the solid matrix is assumed to be linear isotropic elastic, undergoing infinitesimal strains $\epsilon =(gradu+gradTu)/2$, such that

^{f}where $D=(gradvf+gradTvf)/2$ is the rate of deformation tensor, and $\tau f=\phi f\tau $. Finally, we assume that the hydraulic permeability tensor is $k=k\u2009I$ where *k* is constant, and the solute diffusivities are $d\iota =d\iota I$ where $d\iota $ is constant.

for each of the solutes.

### 3.1 Gravitational Sedimentation.

*ι*($\kappa \u0303\iota =1$), subjected to a gravitational specific body force

*g*along the downward positive $z\u2212$direction. As there is no solid matrix, we may set $\phi s=0,\u2009\phi f=1$, $\lambda s=\mu s=0$ (no solid stiffness), $k\u22121=0$ (infinite permeability), and $d\iota =d0\iota I$ (no solute hindrance). Under the conditions of this problem the solvent is at rest ($w=0,\u2009vf=0$), since we assume that the fluid column rests on the ground. Now, the solid momentum balance Eq. (3.3) becomes trivially satisfied ($0=0$), while the solvent and solute momentum Eqs. (3.3)–(3.5) reduce to

*sedimentation coefficient*of solute

*ι*. To complete this set of equations, we also need the mass balance equation for solute

*ι*, reduced from Eq. (2.8) to

*h*is its initial height along

*z*, and the error function is defined by

and the above solution was obtained under the assumption that $Pe\u226b1$. Substituting the solution of Eq. (3.10) for $c\iota $ into the first differential equation of Eq. (3.7) to calculate $jz\iota $ and the resulting expression into Eq. (3.6) can produce a closed form solution for $Jf(z,t)$, though it is too cumbersome to present here. When neglecting the solute flux $jz\iota $, the solution is simply $Jf(z)=exp(\u2212\rho TrfKgz)\u22481\u2212\rho TrfKgz$, such that the effective fluid pressure $p\u0303\u2248\rho Trfgz$ varies linearly with height. The actual fluid pressure is $p(z,t)\u2248p\u0303(z)+R\theta c\iota (z,t)$.

### 3.2 Steady Viscous Electrolyte Flow Over a Charged Porous Layer.

Consider a viscous electrolyte consisting of a solution of two monovalent counterions $\iota =+,\u2212$, flowing above a hybrid multiphasic domain with identical interstitial fluid, whose porous solid matrix has a referential fixed-charge density $crF$ (which we assume to be positive in this illustrative problem). The viscous fluid domain has a height *a* along the positive $y\u2212$direction, while the multiphasic domain has a height *b* along the negative $y\u2212$direction; both two-dimensional domains are infinite along *x*, producing a plane strain analysis. At *y *=* a* the fluid is entrained along *x* at a steady speed *v*_{0}. At $y=\u2212b$ the solid matrix is anchored to a rigid substrate, such that the solid displacement **u** satisfies $u(y=\u2212b)=0$ and the fluid velocity is similarly zero due to the no-slip condition. Due to electroneutrality, the solute concentrations in the viscous electrolyte layer are equal to each other and given by $c+=c\u2212=c0$. The fluid solutions in both domains are assumed to be ideal.

**w**, $j\iota $) are only along

*x*while the dependent variables vary only along the $y\u2212$direction, and we evaluate the $x\u2212$ and $y\u2212$components of those vectorial equations in the absence of external body forces

for the solute momentum balance. The electroneutrality condition requires that $c+\u2212c\u2212+crF=0$. In the fluid layer we set $crF=0$, $\phi f=1,\u2009\lambda s=\mu s=0,\u2009k\u22121=0$ and $d\iota =d0\iota $. In that layer, the solution to the above momentum balance equations, subject to boundary conditions as outlined in Sec. 2.7, produces $Jf=(1\u22122R\theta c0/K)\u22121$ ($p\u0303=\u22122R\theta c0$), $c+=c\u2212=c0$, *p *=* *0, $wx=vxf=(v0\u2212w0)y/a+w0$, and $jx+=jx\u2212=c0wx$, where *w*_{0} is the slip velocity at the interface between the fluid and mixture layers. The electric potential is constant and set to zero (*ψ* = 0).

*y*is produced by the Donnan pressure

*f*increases when the solid matrix hinders solute transport ($d\iota /d0\iota <1$), thus reducing

*ζ*. By satisfying the jump condition (2.107), we find that the slip velocity is given by

*y*=

*0, and the boundary condition*

*u*= 0 at $y=\u2212b$ to produce

_{x}*y*=

*0 the mixture shear traction is $\mu f(v0\u2212w0)/a$ (which is also the shear traction in the fluid layer according to Eq. (2.93)) whereas the solid matrix shear traction is*

This solution shows that the electric potential *ψ* in the charged mixture layer, given in Eq. (3.18), is produced by the fixed-charge density only, regardless of the occurrence of motion. In other words, there are no streaming potentials produced under the conditions of this analysis. However, the electric potential in the charged mixture is affected by the concentration of solutes in the fluid layer flowing above it. Though the above analysis is strictly valid under steady-state only, it remains evident that an evolving or fluctuating concentration *c*_{0} of monovalent counterions in the fluid layer can alter the solute concentrations *c*^{–} and *c*^{+} in Eq. (3.16), and thus the electric potential *ψ* of Eq. (3.18) in the charged mixture. When the charged layer represents the glycocalix of an endothelial cell, it is conceivable that transient alterations in electric potential *ψ* produced by fluctuating salt ion concentrations may activate voltage-gated channels in the associated cell membrane.

## 4 Discussion

The objective of this study was to formulate the theoretical foundations of a framework for multiphysics using a hybrid multiphasic mixture theory. The motivation for such a formulation is that the standard multiphasic mixture frameworks [7,8,11,12,46,53], though versatile and applicable to a wide range of problems in biomechanics and biophysics, incorporate neither dynamics of viscous fluids nor fluid compressibility. We have recently shown that fluid compressibility considerably simplifies the finite element implementation of a computational fluid dynamics solver [41], a fluid-solid structure interaction solver [42], and a fluid-biphasic structure interaction solver [44,45]. In anticipation of incorporating solutes and reactive processes in future finite element implementations, in this study it was necessary to first derive a set of governing equations consistent with the underlying mechanics framework of those earlier formulations.

In the standard multiphasic framework all constituents, including the solvent, are idealized as intrinsically incompressible materials. Consequently the fluid pressure *p* in that standard framework is a Lagrange multiplier which emerges from enforcing this incompressibility constraint in the entropy inequality. It follows that those standard formulations do not include a kinematic state variable to model solvent compressibility. In contrast, as shown here and in our recent hybrid biphasic formulation, the fluid pressure *p* is a function of state. Since we chose to model isothermal processes only, the pressure does not depend on temperature. In this multiphasic framework that incorporates solutes, *p* depends on *J ^{f}*, which is a measure of the fluid volumetric strain, and on solute concentrations, the latter being consistent with concepts of classical physical chemistry.

Indeed, a significant challenge of the present hybrid multiphasic formulation was the selection of state variables for the specific free energies of each constituent, as presented in Eq. (2.49). In a fully general formulation one would normally assume, from the principle of equipresence, that all of the $\psi \alpha $'s should depend on all state variables (unless precluded by the entropy inequality). However, such a general formulation would be too complex to be of practical use. Instead, we had to make careful choices that would maintain the same level of generality as standard multiphasic frameworks as well as classical physical chemistry of dilute solutions, in an isothermal framework that remains suitable for modeling biological tissues and processes, while accommodating the compressibility and viscosity of the solvent.

With these constitutive assumptions, we were able to recover the standard multiphasic split of the fluid pressure *p* into a “hydraulic” contribution $p\u0303$ (which we call *effective* fluid pressure) and an osmotic contribution $R\theta \Phi \u2211\iota c\iota $, as per Eq. (2.84). Here, the effective pressure $p\u0303$ is only a function of *J ^{f}*, even though

*p*depends on

*J*and solute concentrations $c\iota $. Moreover, we recover a logarithmic dependence of the solvent chemical potential

^{f}*μ*on

^{f}*J*when we pick the simplest finite strain constitutive model for $p\u0303(Jf)$, as shown in Eqs. (2.101) and (2.103). This logarithmic dependence is analogous to the starting point for formulating constitutive relations for the chemical potentials of constituents in ideal and nonideal solutions, as outlined in Sec. 2.7.

^{f}In the hybrid biphasic theory, where the mixture consists of a porous solid with intrinsically incompressible skeleton and compressible fluid, we had shown that *J ^{f}* is continuous across interfaces between such domains [44], consistent with the finding from standard biphasic theory that the effective fluid pressure $p\u0303$ is continuous. However, as is well known in physical chemistry, the fluid pressure

*p*of solvent-solute mixtures is not generally continuous across an interface; instead, it is the chemical potential of the solvent and the electrochemical potential of solutes that satisfy this continuity requirement [7,12,69]. The analysis of jump conditions for the hybrid multiphasic theory presented here demonstrates that the same continuity conditions prevail when the solvent is taken to be compressible. As per our recent studies [10,39], we may substitute the continuity of the solvent chemical potential with continuity of $p\u0303$ or, equivalently, continuity of

*J*. The continuity of the electrochemical potential may be substituted with continuity of the effective solute concentrations $c\u0303\iota $, as given in Eq. (2.81), which was also the case in the standard multiphasic framework [10,12,39].

^{f}The hybrid multiphasic framework presented in this study requires the solution of the momentum balance equations for the solid, solvent and solutes in Eqs. (2.85)–(2.87), the mass balance equations for the solid constituents [39], the mass balance for solutes in Eq. (2.8), and the mixture mass balance in Eq. (2.9). The solution of these equations requires additional constitutive models for the porous solid stress $\sigma e$, the viscous fluid stress $\tau $, the hydraulic permeability **k**, the solute diffusivities $d\iota $ and $d0\iota $, the reactive solid mass density supplies $\rho \u0302r\sigma $, the solute molar supplies $c\u0302\iota $, the osmotic coefficient $\Phi $ and the effective solute solubilities $\kappa \u0302\iota $. For most practical problems, these solutions must be obtained using numerical methods.

However, to illustrate the tractability of the equations and boundary conditions presented in this formulation, we provided two analytical solutions for relatively simple canonical problems that may be relevant to important biomechanics problems. The first of these is the one-dimensional sedimentation analysis under the action of gravity, which could be used to model the sedimentation of saliva droplets in air. The reduction of the governing hybrid multiphasic formulation to the conditions of this particular application recovered the classical Smoluchowski convection-diffusion equation, whose mathematical solution summarized above has been reported for various applications before, such as electrophoresis of charged solutes under the action of an electric field. The minor novelty of this presentation is that it included a compressible solvent, though the final solution was obtained assuming nearly incompressible solvent behavior. The second illustrative problem consisted of analyzing planar Couette flow of an electrolyte solution over a charged porous-deformable hybrid multiphasic layer. As far as we know, this solution is novel though it simply represents the superposition of the classical Donnan equilibrium solution [68] with Couette flow over a porous layer [67,70] which is deformable [71]. The solution hints at the possibility that fluctuations in solute concentrations in the flowing electrolyte layer could potentially trigger voltage-gated ion (e.g., calcium) channels in cell membranes due to alterations of the electric potential in the charged glycocalyx.

The primary limitations of this hybrid multiphasic theory are that it is strictly valid only for isothermal processes, and interstitial fluid solutions where solutes occupy a negligible volume fraction. A more subtle limitation is that the specific free energies of all constituents are independent of the (potentially evolving) solid composition $\rho r\sigma $. This constraint, which is inherent to our prior formulations of reactive standard multiphasic frameworks [39,46,47], is not particularly restrictive, though it needs to be mentioned for completeness.

In summary, the formulation of this hybrid multiphasic mixture theory makes it possible to address a wider range of problems that are important in biomechanics and mechanobiology. The next step is to provide computational implementations of this broader framework in open-source software, to enhance the accessibility of this framework across the wider biomechanics and biophysics communities.

## Disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the National Science Foundation.

## Funding Data

Division of Graduate Education, U.S. National Science Foundation (Grant No. NSF GRFP DGE-16-44869; Funder ID: 10.13039/100000082) .

National Institute of General Medical Sciences, U.S. National Institutes of Health (Award No. R01GM083925; Funder ID: 10.13039/100000057).

## References

**143**(9), p.