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 [14], which had led to significant advances in characterizing the mechanics of biological tissues and cells [510]. 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 [1113]. Variations of multiphasic theory have been used in diverse fields such as articular cartilage mechanics [1417], osmotic loading of cells [18], solute transport induced by dynamic loading of porous tissues [1921], 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 [3537]. 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,3840], 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.

We denote generic constituents of a multiphasic mixture by α, with α=s for the solid, α=f for the fluid solvent, and α=ι for the solutes. The definitions of the apparent density ρα, the volume fraction φα, and true density ρTα for each constituent α are given in prior treatments [12,44] and not repeated here. Solutes occupy a negligible volume fraction, φι1 and their true densities ρTι are constant, due to their intrinsic incompressibility. The mixture saturation condition (the absence of voids in the mixture) implies that φs+φf1. 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 Jf relating the solvent current and reference states
Jf=ρTrfρTf=φfρfρTrf=dVfdVrf
(2.1)

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

As described in our standard reactive multiphasic framework [39], the solid matrix 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
ρs=σρσ,φsρTs=σφσρTσ,φs=σφσ
(2.2)

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].

2.2 Motion and Mass Balance.

In mixture theory [1,3,46] each constituent has its own motion, χα(Xα,t), from which the velocity vα may be evaluated [12,44]. As done conventionally, we define the hybrid multiphasic domain on the boundaries of its porous solid constituent. The solid deformation gradient Fs=χs/Xs thus describes the deformation of the porous solid matrix and the solid relative volume Js=detFs is a measure of its volumetric strain. From standard kinematic arguments [12,44], we know that
DsJsDt=Jsdivvs
(2.3)

where DαfDt=ft+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 Jf 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 Jf 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.

When accounting for mass exchanges from reactive processes, the differential statement of mass balance for each constituent is given by
ραt+div(ραvα)=DαραDt+ραdivvα=ρ̂α
(2.4)
where ρ̂α is the mass density supply to α due to reactions with all other constituents. The temporal evolution of the referential solid mass densities ρrσ=Jsρσ and associated solid volume fractions φrσ=Jsφσ 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
φf=1φs=1φrsJs=dVfdV
(2.5)

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

For the solute constituents, it is common to substitute their mass concentration ρι in the mixture with their molar concentration cι in the interstitial fluid [7,39]
cι=ριφfMι,ĉι=ρ̂ιφfMι
(2.6)
where Mι is the molar mass of solute ι, a constant value. Here, we have also provided a similar expression for the molar concentration supply ĉι. In addition, it is also convenient to define the molar solute flux relative to the solid
jι=φfcι(vιvs)
(2.7)
where vι is the solute velocity. The multiplication by φf is needed to account for the fact that this molar flux is normalized by the mixture (not the fluid) area. Using Eqs. (2.4), (2.6), and (2.3), we obtain an expression for the mass balance of the solute constituents
1JsDs(Jsφfcι)Dt+divjι=φfĉι
(2.8)
Substituting the true density into the mass balance (2.4), simplifying them by recalling that ρTσ and ρTι are all constant, taking the sum of the resulting expressions over all α, and making use of the mixture saturation condition αφα=1, produces the mixture mass balance equation
DfJfDt=Jfφf(div(vs+w)αρ̂αρTα)
(2.9)
where
w=φf(vfvs)
(2.10)

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].

Remark 6. In a chemical reaction, the mass density supplies ρ̂α may be related to the molar production rate ζ̂ via ρ̂α=φfMαναζ̂, where να is the net stoichiometric coefficient (product minus reactant) of α in the reaction [39]. Therefore we may rewrite
αρ̂αρTα=φfζ̂αναVα
(2.11)

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).

2.2.1 Electroneutrality Condition.

We assume that electroneutrality is satisfied at all times in the hybrid multiphasic domain, as done in the standard multiphasic theory [10,12], such that
αzαραMα=0,αzαρ̂αMα=0
(2.12)

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.

Assuming that the solvent is electrically neutral (zf=0), and using Eq. (2.6), the electroneutrality condition may be rewritten as
cF+ιzιcι=0
(2.13)
where
cFzsρsφfMs=σzσρσφfMσ
(2.14)
is the fixed charge density of the porous solid matrix, which may evolve with reactions involving the solid constituents. As shown previously [39], we may relate cF to the referential fixed charge density crF via
cF=1φrsJsφrscrF
(2.15)
Multiplying (2.8) by the constant zα/Mα, taking the summation of the resulting equations over all α, and making use of Eqs. (2.6) and (2.12) produces the constraint
0=ιzιdivjι
(2.16)

or equivalently divIe=0, where IeFcιzιjι is the electric current density and Fc 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.

The axiom of linear momentum balance for mixture constituents is given by [12,49]
ραaα=divσα+ραbα+p̂α
(2.17)
Here, σα is the apparent Cauchy stress of α, bα is the body force per mass acting on constituent α, p̂α is the momentum supply per volume to α due to momentum exchanges with other mixture constituents, and aα is the acceleration of α, given by the material time derivative of vα 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
σ=σIαραuαuα
(2.18)
where
σI=ασα
(2.19)
is the inner part of the mixture stress, and
uα=vαv
(2.20)
is the 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
αp̂α+ρ̂αuα=0
(2.21)

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].

2.4 Energy Balance.

Though we plan to use this hybrid multiphasic theory for isothermal processes, we must still formulate the axiom of energy balance for mixture constituents in order to properly formulate the entropy inequality and the interfacial jump conditions needed to prescribe suitable boundary conditions in well-posed hybrid multiphasic problems. The relevant equations are an extension of the equations from our previous hybrid biphasic theory [44]. The axiom of energy balance for each constituent α is repeated here [12,49]
ραDαεαDt=σα:Lαdivqα+ραrα+ε̂α
(2.22)
where εα is the specific internal energy of α, Lα=gradvα is the velocity gradient, qα is the heat flux through α, rα is the specific heat supply to α from sources that are not explicitly modeled, and ε̂α 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
αρ̂α(εα+12uα·uα)+p̂α·uα+ε̂α=0
(2.23)

2.5 Entropy Inequality.

As commonly done for mixture models of biological tissues, we assume that all constituents share the same temperature θ [5,7,44,51]. The mixture entropy inequality is given in the form of the Clausius-Duhem inequality as
αραDαηαDt+divqαθραrαθ0
(2.24)
We may expresss the specific internal energy εα of a constituent as
εα=ψα+θηα
(2.25)
where ψα is the specific free energy, and ηα is the specific entropy of constituent α. The free energy density of α is then defined by Ψα=ραψα where the inner part of the mixture free energy density is given by ΨI=αΨα. The referential free energy density of the mixture is
ΨrI=JsΨI
(2.26)
We may similarly define the entropy density of α as Hα=ραηα and the mixture entropy density as H=αHα, where the referential mixture entropy density is
Hr=JsH
(2.27)
Previous mixture theory presentations [12,46] showed that by using Eqs. (2.22), (2.23), and (2.24)–(2.27), the axiom of entropy inequality for a reactive mixture may take the form
1Js(DsΨrIDt+HrDsθDt)+(σs+αs(ΨαzαMαραFcψ)I):Ls+αs(σα+(zαMαραFcψΨα)I):Lα(p̂sαs(gradΨα+HαgradθzαMαFcψgradρα))·usαs(p̂α+gradΨα+HαgradθzαMαFcψgradρα)·uα1θgradθ·αqαα12ρ̂αuα·uα0
(2.28)
This expression also incorporates the electroneutrality constraint (2.16), rewritten in the form
αzαMαdiv(ραuα)=αzαMα(ραLα:I+gradρα·uα)=0
(2.29)

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].

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 Jf 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.

Using the chain rule of differentiation on ΨrI(Fs,Jf,Lf,ρrσ,ρrι,uα) produces
DsΨrIDt=ΨrIFs:DsFsDt+ΨrIJfDsJfDt+σΨrIρrσDsρrσDt+ιΨrIρrιDsρrιDt+ΨrILf:DsLfDt+ΨrIus·DsusDt+ΨrIuf·DsufDt+ιΨrIuι·DsuιDt
(2.30)
In this expression we may substitute ρ̂rσ=Jsρ̂σ for Dsρrσ/Dt based on the mass balance for the solid, as well as the standard kinematic identity
DsFsDt=Ls·Fs
(2.31)
expressed here in the spatial frame. Using the definition of Jf 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
DsJfDt=JfφsφfI:Ls+JfI:Lf+Jfρfgradρf·(ufus)Jfφfσρ̂σρTσ
(2.32)
and substituted into Eq. (2.30). Similarly, the solute mass balance (2.8) may be rearranged in the form
DsρrιDt=ρ̂rιJsgradρι·(uιus)ρrιI:(LιLs)
(2.33)
to be substituted into Eq. (2.30). The resulting expression in Eq. (2.30) may be substituted into the entropy inequality (2.28) and, upon regrouping terms, produces
1Js(HrDsθDt+ΨrILf:DsLfDt+αΨrIuα·DsuαDt)+(σs1JsΨrIFs·(Fs)T+ρTfμfI+αs(Ψαραμ̃α)I):Ls+(σf(Ψfρfμf)I):Lf+ι(σι(Ψιριμ̃ι)I):Lι(p̂sαs(gradΨα+Hαgradθμ̃αgradρα))·usαs(p̂α+gradΨα+Hαgradθμ̃αgradρα)·uα1θgradθ·αqααρ̂α(μα+12uα·uα)ρTfμfσρ̂σρTσ0
(2.34)
In this expression note that
μαΨrIρrα
(2.35)
is the chemical potential of constituent α, while its electrochemical potential is given by
μ̃αμα+zαMαFcψ
(2.36)
When constituent α is electrically neutral, zα=0 and μ̃α=μα. Using the definition of Eq. (2.35) for α=f, along with Eq. (2.1) and ρrf=Jsρf, we also note that the chemical potential of the fluid solvent is related to ΨrI/Jf via
μfΨrIρrf=JfJsρfΨrIJf
(2.37)

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

For arbitrary processes (i.e., arbitrary variations in the variables of state Dsθ/Dt,gradθ,Lα, uα,DsLf/Dt,Dsuα/Dt), the inequality in Eq. (2.34) is satisfied unconditionally if and only if the following constraints are enforced
Hr=0
(2.38)
ΨrILf=ΨrIuα=0
(2.39)
and
αqα=0
(2.40)
We also obtain
σs=1JsΨrIFs·(Fs)TρTfμfIαs(Ψαραμ̃α)I
(2.41)
and
σι=(Ψιριμ̃ι)I
(2.42)
leaving the residual dissipation statement
τf:Lfαp̂dα·uααρ̂α(μα+12uα·uα)ρTfμfσρ̂σρTσ0
(2.43)
Here
τfσf(Ψfρfμf)I
(2.44)
represents the dissipative (viscous) part of the fluid stress, whereas
p̂ds=p̂sαs(gradΨα+Hαgradθμ̃αgradρα)
(2.45)
and
p̂dα=p̂α+gradΨα+Hαgradθμ̃αgradρα,αs
(2.46)
represent the dissipative parts of the momentum supplies. From Eq. (2.38) we recover the standard finding that the mixture entropy is zero in an isothermal framework. We can easily verify from these relations and Eq. (2.21) that the dissipative momentum supplies also satisfy αp̂dα=αp̂α, such that Eq. (2.17) now produces
αp̂dα+ρ̂αuα=0
(2.47)
Substituting Eqs. (2.41), (2.42), and (2.44) into Eq. (2.19) produces
σI=ρTfμfI+1JsΨrIFs·(Fs)T+τf
(2.48)

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.

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, ψβ=ψβ(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, Js 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ι.

Since the solvent is compressible, we also assume that ψf depends on Jf. Similarly, since the porous solid is deformable, we assume that ψσ depends on Fs. Finally, we adopt the simplifying assumption that ψσ does not depend on any of the (potentially evolving) solid concentrations ρrσ, implying that solid composition influences the referential free energy density of each solid constituent σ only based on the proportionality relation Ψrσ=ρrσψσ. With these constitutive assumptions, it follows from Eq. (2.26) that
ΨrI(Fs,Jf,ρrσ,ρrι)=σρrσψσ(Fs,cι)+ρrfψf(Jf,cι)+βρrβψβ(cι)
(2.49)
where ρrf=Jsρf and Jf depend on each other as per Eq. (2.1). Accordingly, using Eq. (2.35) we get μσ=ψσ for solid constituents. Similarly, using Eq. (2.37), we find that
μf=ψfJfψfJf
(2.50)
for the fluid chemical potential and
σI=ρTrfψfJfI+1JsΨrsFs·(Fs)T+τf
(2.51)
for the mixture stress. We may now define the pressure in the compressible fluid as
p(Jf,cι)ρTrfψfJf
(2.52)
which is similar to the expression for fluid pressure in our compressible isothermal fluid [41] and hybrid biphasic formulations [44,45]. Combining Eq. (2.50) with Eq. (2.52) produces an alternative form for μf
μf=ψf+pρTf
(2.53)

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.

Substituting the constitutive model Eq. (2.49) of Eq. (2.26) into Eq. (2.35), we also find that the chemical potential of solute β takes the form
μβ=ψβ+Ψrsρrβ+ρrfψfρrβ+ιρrιψιρrι
(2.54)
In the expression for the inner part of the mixture stress, Eq. (2.51), the elastic stress in the porous solid may be defined as
σe1JsΨrsFs·FsT
(2.55)
In both expressions above, the solid referential free energy density is given by
Ψrs=σρrσψσ
(2.56)
Thus, reactions involving the solid constituents may produce an evolution in Ψrs due to the evolving solid compositions ρrσ. From these expressions, we can simplify (2.41), (2.44), and (2.42) to get
σs=φspI+σeι(Ψιριμ̃ι)I
(2.57)
and
σf=φfpI+τf
(2.58)
while σι does not simplify further. Now, the inner part of the mixture stress in Eq. (2.19) may be rewritten as
σI=ασα=pI+σe+τf
(2.59)
Substituting these stresses and momentum supplies into each of the constituent momentum balances (2.17) and assuming isothermal conditions (gradθ = 0) produces
ρsas=gradp+ρfgradμf+divσe+ρsbs+p̂ds+ιριgradμ̃ι
(2.60)
for the solid
ρfaf=ρfgradμf+divτf+ρfbf+p̂df
(2.61)
for the solvent, and
0=ριgradμ̃ι+ριbι+p̂dι
(2.62)

for the solutes, where inertia terms have been neglected.

2.6 Constitutive Relations.

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:Lf0 and αp̂dα·uα0, which are sufficient (but not necessary) conditions to satisfy the residual dissipation statement (2.43).

2.6.1 Dissipative Momentum Supplies.

Following our previous presentation [12], for a reactive mixture we can let
p̂dα=ραρβρ̂βuβ+βαfαβ(uβuα)
(2.63)
where fαβ=fαβ(Fs,Jf,Lf,ρrσ,ρrι,uα) is the diffusive drag tensor. To satisfy the constraint (2.47) for arbitrary diffusive velocities, we must enforce
fαβ=fβα
(2.64)

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 τf0.

In this hybrid multiphasic mixture model, we may account for frictional interactions between the solvent and solid constituents as well as solvent viscosity, to reproduce Darcy-Brinkman type of friction [52]. We may also model frictional interactions between the solvent and solutes and the solid and solutes [53], such as those that result in Fickian diffusion [54] and hindrance mechanisms for solute transport in porous media [55]. However, frictional interactions between solutes of the same species (solute viscosity) and between solutes of different species are neglected, since the interstitial fluid solution is dilute. Since reactions are assumed to maintain dilute conditions, we may also neglect the reactive term on the right-hand side of Eq. (2.63) relative to the last term. As shown previously [7,12,53], these simplifying assumptions produce relations between the diffusive drag tensors fαβ and the hydraulic permeability tensor k of the solvent through the porous solid matrix, the solute diffusivity d0ι in free solution, and the (hindered) solute diffusivity tensor dι in the porous-permeable mixture according to
ffs=(φf)2k1
(2.65)
fιf=Rθφfcι(d0ι)1I
(2.66)
and
fιs+fιf=Rθφfcι(dι)1
(2.67)

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) [5659], 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.

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
ψι=ψ0ι+RθMιlnxι
(2.68)
where xι is the mole fraction of solute ι in the solution. Similarly, the free energy of the solvent is given by
ψf=ψ0f+RθMflnxf
(2.69)
where xf is the mole fraction of the solvent. In a solution, mole fractions satisfy αsxα=1. Thus, the mole fraction of the solvent may also be evaluated as xf=1ιxι, 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ψ0fRθMfιxι
(2.70)
where ψ0f is the value of ψf when the solution contains no solutes. Thus, ψι depends only on its own xι, whereas ψf depends on all xι's. Note that xf is not an independent state variable in this context. For a dilute solution, we can express the mole fraction of solute ι in terms of its molar concentration cι as
xιMfρTfcι
(2.71)
Thus, the solute specific free energy in Eq. (2.68) may be rewritten as
ψι(cι)=ψ0ι+RθMιlnMfc0ρTfμ0ι+RθMιlncιc0μ0ι+RθMιlncιc0
(2.72)
where c0 is an arbitrary constant defining the solute standard state (typically, c0=1M) and μ0ι is a constant that represents the value of ψι when cι=c0. Similarly, the solvent specific free energy in Eq. (2.70) may be rewritten as
ψfψ0fRθρTfιcι
(2.73)
Using the expression for μι in Eq. (2.35), along with Eq. (2.49) (where ρrσ=0 for this fluid mixture), we get
μι=ΨrIρrι=ΨIρι=ρι(ρfψf+αfραψα)=ρfψfρι+ψι+ριψιρι
(2.74)
Now, we use (2.72)(2.73), ρf=ρTf=Mfcf and ρι=Mιcι to reduce this expression to
μι=μ0ι+RθMιlncιc0
(2.75)
Similarly, using the expression for μf in Eq. (2.53), we get
μf=ψ0f+1ρTf(pRθιcι)
(2.76)

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.
For nonideal solutions and mixtures that include solid constituents, we follow a standard approach in physical chemistry that generalizes the formulation from ideal solutions such that
μι=μ0ι+RθMιlnaι
(2.77)
is the generalization of Eq. (2.72), where aι is the activity of solute ι and μ0ι is the chemical potential in the standard state when aι=1. For nonideal solutions, we can express aι as
aι=γιcικιc0
(2.78)
where the activity coefficientγι(Fs,Jf,ρrσ,cβ) 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), κι 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κι1. For ideal mixtures, we let γι=1 and for perfect solubility we let κι=1. Since γι and κι always appear as a ratio in these equations, we may define the effective solubility [61] as
κ̂ικιγι
(2.79)
Now, the electrochemical potential of solute ι may be written as
μ̃ι=μ0ι+RθMιlnc̃ιc0ι
(2.80)
Here, the effective concentration is defined as
c̃ιcικ̃ι
(2.81)
where
κ̃ι=κ̂ιexp(zιFcψRθ)
(2.82)

is the partition coefficient of ι relative to an ideal solution [62,63].

For the solvent in a mixture, we generalize the expression of Eq. (2.76) to
μf=ψ0f+1ρTf(pRθΦιcι)
(2.83)
Here, Φ=Φ(Fs,Jf,ρrσ,cι) is the 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
p̃=pRθΦιcι
(2.84)

whereas RθΦιcι represents the “chemical” or “osmotic” contribution to the fluid pressure.

2.6.3 Implication for Momentum Balances.

Finally, taking into account constitutive relations from Secs. 2.6.1 and 2.6.2, we can rewrite the momentum balances (2.60), (2.61), and (2.62) as
ρsas=gradp+ρfgradμf+ιριgradμ̃ι+divσe+ρsbs+φfk1·w+Rθι((dι)1(d0ι)1I)·jι
(2.85)
ρfaf=ρfgradμf+divτf+ρfbf(φfk1+Rθιcι(d0ι)1I)·w+Rθι(d0ι)1I·jι
(2.86)
and
0=ριgradμ̃ι+ριbιRθ(dι)1·jι+Rθcι(d0ι)1w
(2.87)

2.7 Jump Conditions Across Reactive Interfaces.

To solve the above governing equations, we must determine the correct set of boundary conditions across the mixture interfaces. The boundary conditions are obtained from the jump conditions on the axioms of mass, momentum and energy balance across an immaterial interface Γ. Various authors have previously reported jump conditions for general mixtures [12,65,66]. The derivations of these jump conditions are similar to our previous derivations for the hybrid biphasic mixture [44], except for the addition of solute constituents. For the axiom of mass balance, the jump condition for each constituent α takes the form
[[ραuΓα]]·n=ρ¯α
(2.88)

where uΓα=vαvΓ, vΓ is the velocity of interface Γ, n is a unit normal on Γ, and ρ¯α is the mass flux of α produced by reactions on Γ. In this notation, for any arbitrary vector v, we have [[v]]·n=v+·n++v·n=(v+v)·n, 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=n. In multiphasic mixtures, we may define the interface Γ on a porous solid boundary, where vΓ=vs and uΓs=0.

For each solute ι, substituting Eqs. (2.6) and (2.7) into Eq. (2.88) produces
[[jι]]·n=c¯ι
(2.89)

where c¯ι=ρ¯ι/Mι 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¯ι>0) or consume (c¯ι<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.

As shown previously for the hybrid biphasic mixture [44], for the nonreactive solvent the jump condition (2.88) reduces to
[[ρTrfJfw]]·n=0
(2.90)
which implies that the component of the solvent mass flux relative to the solid is continuous in the direction normal to Γ. We similarly formulate a condition for the tangential viscous fluid velocity on Γ, by following the approach of Hou et al. [67], who substituted n in Eq. (2.90) with the unit vector s tangent to Γ [44], producing
[[ρTrfJfw]]·s=0
(2.91)
Hou et al. [67] called this condition the “pseudo-no-slip” condition, since it reduces to the standard no-slip condition when one side is a viscous fluid (ρf0) and the other is an impermeable solid (ρf=0). Combining normal and tangential components in Eqs. (2.90) and (2.91) produces
[[ρTrfJfw]]=0
(2.92)
The jump condition from the axiom of linear momentum balance can be simplified similarly to that of the hybrid biphasic mixture [44], where
[[σI]]·n=[[pI+τf+σe]]·n0
(2.93)

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

The jump condition resulting from the axiom of energy balance is also reduced according to our previous formulation [44], but here we also include the solute stress (2.42) to obtain
n·[[τfρf·ρTrfJfw]]+[[ρTrfJfμfw]]·n+ι[[ριμ̃ιuΓι]]·n=0
(2.94)
Based on the continuity of the solvent mass flux in Eq. (2.90), we can take ρTrfw/Jf out of the double brackets in this expression. The same can be done for ριuΓι for the solutes if solute ι is not involved in a reaction on Γ (thus, ρ¯ι=0 in Eq. (2.88)). In that case, Eq. (2.94) may be simplified to
n·[[τfρf·ρTrfJfw]]+[[μf]]ρTrfJfw·n+ι[[μ̃ι]]ριuΓι·n=0
(2.95)
This jump condition may be unconditionally satisfied on a permeable interface if the following sufficient (but not necessary) jump conditions are satisfied:
[[τfρf]]·n=0
(2.96)
which requires a specific measure of the fluid viscous traction to be continuous across Γ
[[μf]]=0
(2.97)
where the fluid mechanochemical potential must be continuous across Γ, and
[[μ̃ι]]=0,nonreactivesolutes
(2.98)

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¯ι as a flux (Neumann) boundary condition on ι, instead of satisfying (2.98) as a Dirichlet boundary condition.

Remark 13. For an ideal solution, as reviewed in Sec. 2.6.2, the solvent specific free energy in Eq. (2.73) may be rewritten as
ψf(Jf,cι)=ψ0f(Jf)RθρTrfJfιcι
(2.99)
where ψ0f is that part of ψf that does not depend on solute concentrations. Substituting this expression into the general relation of Eq. (2.50) produces
μf=ψ0fJfdψ0fdJf
(2.100)

This result demonstrates that μf is exclusively a function of Jf when the mixture is an ideal solution. It follows from Eq. (2.76) that the effective pressure p̃=pRθιcι of an ideal solution is also only a function of Jf, even though the actual fluid pressure p depends on Jf and the solute concentrations, as per (2.52).

For nonideal mixtures, combining (2.83) and (2.84) produces μf=ψ0f+p̃/ρTf, where ψ0f and ρTf are only functions of Jf. Based on the above remark regarding ideal solutions, and the fact that p̃ is the nonosmotic part of the pressure p as given in Eq. (2.84), it follows that p̃(Jf) for nonideal mixtures must also depend only on the fluid volume ratio, implying that μf=μ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 Jf
p̃(Jf)=ρTrfdψ0fdJf=K(1Jf1)
(2.101)
implying that
ψ0f(Jf)=KρTrf(Jf1lnJf)
(2.102)
and
μf=KρTrflnJf=KρTrflnρTfρTrf
(2.103)
Thus, using the constitutive model of Eq. (2.101), the solvent chemical potential exhibits the familiar logarithmic form also employed for solutes as shown above. As a consequence of this exclusive dependence on Jf, the continuity of the solvent chemical potential μf across Γ, expressed in Eq. (2.97), implies continuity of Jf, on the understanding that the same solvent species is present on both sides of Γ, thus
[[Jf]]=0
(2.104)
This condition implies that the effective pressure p̃ must also be continuous across Γ
[[p̃]]=0
(2.105)
as also found in the standard multiphasic formulation [12]. Using Eq. (2.104), we can simplify Eqs. (2.92) and (2.96) to, respectively, obtain
[[w]]=0
(2.106)
implying that the volumetric flux of the viscous fluid is continuous across Γ, and
[[τfφf]]·n=0
(2.107)
As with the standard multiphasic formulation [12], we can also simplify (2.98) to produce continuity of the effective solute concentration across Γ
[[c̃ι]]=0
(2.108)

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 τf=φfτ where τ would be the viscous stress in the free fluid (absent a porous matrix, when φf=1). Thus, we could use standard constitutive relations for the viscous stress τ of Newtonian or non-Newtonian fluids and adapt those models to a multiphasic mixture where τf is evaluated from φfτ.

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.

We use the governing equations in Eqs. (2.85)–(2.87), specialized to the conditions of the two models analyzed below. For quasi-static or steady-state analyses we let as=af=0. We use ideal physicochemical constitutive relations, including (2.80) for the solute electrochemical potential where we substitute (2.81) and (2.82) with κ̂ι=1. We also use (2.83) for the solvent chemical potential, where we set Φ = 1. We use the relation of Eq. (2.84) to relate p̃ and p, where we insert the constitutive model of Eq. (2.101) for p̃, thus reducing the solvent chemical potential to the form given in Eq. (2.103). We employ (2.1) to solve for ρf and Eq. (2.6) to solve for ρι. When present, the solid matrix is assumed to be linear isotropic elastic, undergoing infinitesimal strains ε=(gradu+gradTu)/2, such that
σe=λs(trε)I+2μsε
(3.1)
As a result, the solid and fluid volume fractions φs and φf, and the fixed charge density remain essentially constant in the porous layer. The viscous fluid in both domains is assumed to be Newtonian, such that
τ=λf(trD)I+2μfD
(3.2)

where D=(gradvf+gradTvf)/2 is the rate of deformation tensor, and τf=φfτ. Finally, we assume that the hydraulic permeability tensor is k=kI where k is constant, and the solute diffusivities are dι=dιI where dι is constant.

Upon substituting all these relations into Eqs. (2.85)–(2.87) we get
0=K(Jf)2φsgradJf+Rθι(φfκ̃ιgradcικ̃ιgradcι)+(λs+μs)grad(divu)+μslapu+φsρTsbs+φfk1w+Rθι((dι)1(d0ι)1)jι
(3.3)
for the solid momentum
0=K(Jf)2φfgradJf+φfJfρTrfbf+φf(λf+μf)grad(divvf)+φfμflapvfφfk1w+Rθι(d0ι)1(jιcιw)
(3.4)
for the fluid momentum, and
jι=φfκ̃ιdιgradcικ̃ι+φfMιdιRθcιbι+dιd0ιcιw
(3.5)

for each of the solutes.

3.1 Gravitational Sedimentation.

We analyze this problem in one dimension of space (the zdirection), for a column of fluid containing a single, electrically neutral solute ι (κ̃ι=1), subjected to a gravitational specific body force g along the downward positive zdirection. As there is no solid matrix, we may set φs=0,φf=1, λs=μs=0 (no solid stiffness), k1=0 (infinite permeability), and dι=d0ιI (no solute hindrance). Under the conditions of this problem the solvent is at rest (w=0,vf=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
0=KJfJfz+ρTrfg+RθJfd0ιjzι
(3.6)
and
jzι=d0ιcιz+sιcιg
(3.7)
where sιMιd0ι/Rθ is known as the 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
cιt+jzιz=0
(3.8)
Substituting the expression for jzι in Eq. (3.7) into the above mass balance produces the standard Smoluchowski partial differential equation that combines diffusion and convection of the solute
cιtd0ι2cιz2+sιgcιz=0
(3.9)
This equation may be solved by assuming that diffusion takes place more slowly than sedimentation, over the time frame of the analysis, implying that |sιgcι||d0ιdcι/dz| on the right-hand-side of the expression for jzι in Eq. (3.7). To get a closed-form solution we also assume that the zdomain is infinite, allowing us to perform a change of variable ζ(z,t)=zsιgt to examine the diffusion of a convecting bolus of solute centered on the bolus. Using the method of Fourier transforms, the final solution for cι(z,t) takes the form
cι(z,t)=c0ι2(erfh(zsιgt)2d0ιt+erfzsιgt2d0ιt)
(3.10)
where c0ι is the initial bolus concentration, h is its initial height along z, and the error function is defined by
erf(x)=2π0xeu2du
(3.11)
A representative plot of this solution at different time points is given in Fig. 1. For this sedimentation-diffusion problem, the sedimentation velocity is sιg whereas the characteristic diffusion velocity is d0ι/h. Therefore, the Peclet number is given by
Pe=hsιgd0ι
(3.12)

and the above solution was obtained under the assumption that Pe1. Substituting the solution of Eq. (3.10) for cι into the first differential equation of Eq. (3.7) to calculate jzι 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ι, the solution is simply Jf(z)=exp(ρTrfKgz)1ρTrfKgz, such that the effective fluid pressure p̃ρTrfgz varies linearly with height. The actual fluid pressure is p(z,t)p̃(z)+Rθcι(z,t).

Fig. 1
Sedimentation problem: normalized solute concentration profile cι/c0ι for Pe=103 versus normalized vertical position z/h, shown at four consecutive values of the nondimensional time t/τ, where τ=h2/d0ι
Fig. 1
Sedimentation problem: normalized solute concentration profile cι/c0ι for Pe=103 versus normalized vertical position z/h, shown at four consecutive values of the nondimensional time t/τ, where τ=h2/d0ι
Close modal

3.2 Steady Viscous Electrolyte Flow Over a Charged Porous Layer.

Consider a viscous electrolyte consisting of a solution of two monovalent counterions ι=+,, 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 ydirection, while the multiphasic domain has a height b along the negative ydirection; 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 v0. At y=b the solid matrix is anchored to a rigid substrate, such that the solid displacement u satisfies u(y=b)=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=c0. The fluid solutions in both domains are assumed to be ideal.

To solve Eqs. (3.3)–(3.5) we assume that the flows (w, jι) are only along x while the dependent variables vary only along the ydirection, and we evaluate the x and ycomponents of those vectorial equations in the absence of external body forces
0=μsd2uxdy2+(φfk+Rθι1d0ι(1dιd0ι)cι)wx0=(λs+2μs)d2uydy2Rθιdcιdy
(3.13)
for the solid momentum balance
0=μfd2wxdy2(φfk+Rθι1d0ι(1dιd0ι)cι)wx0=dJfdy
(3.14)
for the fluid momentum balance, and
jxι=dιd0ιcιwx0=ddy(cικ̃ι)
(3.15)

for the solute momentum balance. The electroneutrality condition requires that c+c+crF=0. In the fluid layer we set crF=0, φf=1,λs=μs=0,k1=0 and dι=d0ι. In that layer, the solution to the above momentum balance equations, subject to boundary conditions as outlined in Sec. 2.7, produces Jf=(12Rθc0/K)1 (p̃=2Rθc0), c+=c=c0, p =0, wx=vxf=(v0w0)y/a+w0, and jx+=jx=c0wx, where w0 is the slip velocity at the interface between the fluid and mixture layers. The electric potential is constant and set to zero (ψ = 0).

In the mixture domain, we also find that Jf=(12Rθc0/K)1 (p̃=2Rθc0); there, by satisfying the jump conditions (2.105) and (2.108), the solute concentrations are given by the Donnan equilibrium solution [68]
c=crF+(crF)2+(2c0)22c+=crF+(crF)2+(2c0)22
(3.16)
the fluid (Donnan) pressure is
p=Rθ((crF)2+(2c0)22c0),
(3.17)
and the electric (Donnan) potential is similarly
ψ=Rθ2Fclncc+
(3.18)
These solutions are homogeneous throughout the domain by0. The solid matrix deformation along y is produced by the Donnan pressure
uy(y)=pλs+2μs(y+b)by0
(3.19)
The horizontal component of the fluid flux is given by
wx(y)=φfvxf=w0sinh[1ζ(yb+1)]sinh1ζby0
(3.20)
where the nondimensional parameter ζ=μf/fb2 controls the thickness of the boundary layer in the vicinity of y0 (Fig. 2). Here
f=φfk+Rθι=+,1d0ι(1dιd0ι)cι
(3.21)
which shows that f increases when the solid matrix hinders solute transport (dι/d0ι<1), thus reducing ζ. By satisfying the jump condition (2.107), we find that the slip velocity is given by
w0=baζtanh1ζ1+φfζbatanh1ζφfv0
(3.22)
Fig. 2
Normalized interstitial fluid flux wx(y)/w0 in the mixture layer (−b≤y≤0) for steady viscous electrolyte flow over a charged porous layer. The nondimensional parameter ζ controls the thickness of the boundary layer near y≲0.
Fig. 2
Normalized interstitial fluid flux wx(y)/w0 in the mixture layer (−b≤y≤0) for steady viscous electrolyte flow over a charged porous layer. The nondimensional parameter ζ controls the thickness of the boundary layer near y≲0.
Close modal
The solute fluxes are given by the first relation in Eq. (3.15), and the solid displacement in the direction of the flow is obtained by solving the first equation in Eq. (3.13), making use of the jump condition (2.93) on the mixture traction at y =0, and the boundary condition ux = 0 at y=b to produce
ux(y)=μfμs(φsba(v0w0)+w0ζtanh1ζ)(yb+1)μfμssinh[1ζ(yb+1)]sinh1ζw0,by0
(3.23)
At the interface y =0 the mixture shear traction is μf(v0w0)/a (which is also the shear traction in the fluid layer according to Eq. (2.93)) whereas the solid matrix shear traction is
σxye|y=0=μsduxdy|y=0=φsμfv0w0a
(3.24)

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 c0 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 Jf, 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 ψα'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̃ (which we call effective fluid pressure) and an osmotic contribution RθΦιcι, as per Eq. (2.84). Here, the effective pressure p̃ is only a function of Jf, even though p depends on Jf and solute concentrations cι. Moreover, we recover a logarithmic dependence of the solvent chemical potential μf on Jf when we pick the simplest finite strain constitutive model for p̃(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.

In the hybrid biphasic theory, where the mixture consists of a porous solid with intrinsically incompressible skeleton and compressible fluid, we had shown that Jf is continuous across interfaces between such domains [44], consistent with the finding from standard biphasic theory that the effective fluid pressure p̃ 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̃ or, equivalently, continuity of Jf. The continuity of the electrochemical potential may be substituted with continuity of the effective solute concentrations c̃ι, as given in Eq. (2.81), which was also the case in the standard multiphasic framework [10,12,39].

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 σe, the viscous fluid stress τ, the hydraulic permeability k, the solute diffusivities dι and d0ι, the reactive solid mass density supplies ρ̂rσ, the solute molar supplies ĉι, the osmotic coefficient Φ and the effective solute solubilities κ̂ι. 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 ρrσ. 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

1.
Truesdell
,
C.
, and
Toupin
,
R.
,
1960
, “
Encyclopedia of Physics
,”
Chapter the Classical Field Theories
, Vol.
III/1
,
Springer-Verlag
,
Berlin
.
2.
Green
,
A. E.
, and
Naghdi
,
P. M.
,
1969
, “
On Basic Equations for Mixtures
,”
Q. J. Mech. Appl. Math.
,
22
(
4
), pp.
427
438
.10.1093/qjmam/22.4.427
3.
Bowen
,
R. M.
,
1976
,
Continuum Physics, Chap. Theory of Mixtures
,
Academic Press
, New York, pp.
1
127
.
4.
Bedford
,
A.
, and
Drumheller
,
D. S.
,
1983
, “
Theories of Immiscible and Structured Mixtures
,”
Int. J. Eng. Sci.
,
21
(
8
), pp.
863
960
.10.1016/0020-7225(83)90071-X
5.
Mow
,
V. C.
,
Kuei
,
S. C.
,
Lai
,
W. M.
, and
Armstrong
,
C. G.
,
1980
, “
Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments
,”
ASME J. Biomech. Eng.
,
102
(
1
), pp.
73
84
.10.1115/1.3138202
6.
Oomens
,
C. W.
,
van Campen
,
D. H.
, and
Grootenboer
,
H. J.
,
1987
, “
A Mixture Approach to the Mechanics of Skin
,”
J. Biomech.
,
20
(
9
), pp.
877
885
.10.1016/0021-9290(87)90147-3
7.
Lai
,
W. M.
,
Hou
,
J. S.
, and
Mow
,
V. C.
,
1991
, “
A Triphasic Theory for the Swelling and Deformation Behaviors of Articular Cartilage
,”
ASME J. Biomech. Eng.
,
113
(
3
), pp.
245
258
.10.1115/1.2894880
8.
Huyghe
,
J. M.
, and
Janssen
,
J.
,
1997
, “
Quadriphasic Mechanics of Swelling Incompressible Porous Media
,”
Int. J. Eng. Sci.
,
35
(
8
), pp.
793
802
.10.1016/S0020-7225(96)00119-X
9.
Humphrey
,
J. D.
, and
Rajagopal
,
K. R.
,
2003
, “
A Constrained Mixture Model for Arterial Adaptations to a Sustained Step Change in Blood Flow
,”
Biomech. Model. Mechanobiol.
,
2
(
2
), pp.
109
–1
26
.10.1007/s10237-003-0033-4
10.
Ateshian
,
G. A.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2013
, “
Multiphasic Finite Element Framework for Modeling Hydrated Mixtures With Multiple Neutral and Charged Solutes
,”
ASME J. Biomech. Eng.
,
135
(
11
), p. 111001.10.1115/1.4024823
11.
Gu
,
W. Y.
,
Lai
,
W. M.
, and
Mow
,
V. C.
,
1998
, “
A Mixture Theory for Charged-Hydrated Soft Tissues Containing Multi-Electrolytes: Passive Transport and Swelling Behaviors
,”
ASME J. Biomech. Eng.
,
120
(
2
), pp.
169
80
.10.1115/1.2798299
12.
Ateshian
,
G. A.
,
2007
, “
On the Theory of Reactive Mixtures for Modeling Biological Growth
,”
Biomech. Model. Mechanobiol.
,
6
(
6
), pp.
423
445
.10.1007/s10237-006-0070-x
13.
Ateshian
,
G. A.
,
2016
, “
Mixture Theory for Modeling Biological Tissues: Illustrations From Articular Cartilage
,”
Studies in Mechanobiology, Tissue Engineering and Biomaterials
,
Springer International Publishing
, Berlin, pp.
1
51
.
14.
Wang
,
C. C.-B.
,
Guo
,
X. E.
,
Sun
,
D.
,
Mow
,
V. C.
,
Ateshian
,
G. A.
, and
Hung
,
C. T.
,
2002
, “
The Functional Environment of Chondrocytes Within Cartilage Subjected to Compressive Loading: A Theoretical and Experimental Approach
,”
Biorheology
,
39
(1–2), pp.
11
25
.https://pubmed.ncbi.nlm.nih.gov/12082263/
15.
Ateshian
,
G. A.
,
Soltz
,
M. A.
,
Mauck
,
R. L.
,
Basalo
,
I. M.
,
Hung
,
C. T.
, and
Lai
,
W. M.
,
2003
, “
The Role of Osmotic Pressure and Tension-Compression Nonlinearity in the Frictional Response of Articular Cartilage
,”
Transp. Porous Media
,
50
(
1/2
), pp.
5
33
.10.1023/A:1020618514874
16.
Ateshian
,
G. A.
,
Chahine
,
N. O.
,
Basalo
,
I. M.
, and
Hung
,
C. T.
,
2004
, “
The Correspondence Between Equilibrium Biphasic and Triphasic Material Properties in Mixture Models of Articular Cartilage
,”
J. Biomech.,
37
(
3
), pp.
391
400
.10.1016/s0021-9290(03)00252-5
17.
Lu
,
X. L.
,
Wan
,
L. Q.
,
Guo
,
X. E.
, and
Mow
,
V. C.
,
2010
, “
A Linearized Formulation of Triphasic Mixture Theory for Articular Cartilage, and Its Application to Indentation Analysis
,”
J. Biomech.,
43
(
4
), pp.
673
679
.10.1016/j.jbiomech.2009.10.026
18.
Albro
,
M. B.
,
Chahine
,
N. O.
,
Caligaris
,
M.
,
Wei
,
V. I.
,
Likhitpanichkul
,
M.
,
Ng
,
K. W.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2007
, “
Osmotic Loading of Spherical Gels: A Biomimetic Study of Hindered Transport in the Cell Protoplasm
,”
ASME J. Biomech. Eng.
,
129
(
4
), pp.
503
510
.10.1115/1.2746371
19.
Albro
,
M. B.
,
Chahine
,
N. O.
,
Li
,
R.
,
Yeager
,
K.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2008
, “
Dynamic Loading of Deformable Porous Media Can Induce Active Solute Transport
,”
J. Biomech.
,
41
(
15
), pp.
3152
3157
.10.1016/j.jbiomech.2008.08.023
20.
Chahine
,
N. O.
,
Albro
,
M. B.
,
Lima
,
E. G.
,
Wei
,
V. I.
,
Dubois
,
C. R.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2009
, “
Effect of Dynamic Loading on the Transport of Solutes Into Agarose Hydrogels
,”
Biophys. J.
,
97
(
4
), pp.
968
975
.10.1016/j.bpj.2009.05.047
21.
Albro
,
M. B.
,
Li
,
R.
,
Banerjee
,
R. E.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2010
, “
Validation of Theoretical Framework Explaining Active Solute Uptake in Dynamically Loaded Porous Media
,”
J. Biomech.
,
43
(
12
), pp.
2267
2273
.10.1016/j.jbiomech.2010.04.041
22.
Chen
,
X.
, and
Sarntinoranont
,
M.
,
2007
, “
Biphasic Finite Element Model of Solute Transport for Direct Infusion Into Nervous Tissue
,”
Ann. Biomed. Eng.
,
35
(
12
), pp.
2145
2158
.10.1007/s10439-007-9371-1
23.
Smith
,
J. H.
, and
García
,
J. J.
,
2011
, “
A Nonlinear Biphasic Model of Flow-Controlled Infusions in Brain: Mass Transport Analyses
,”
J. Biomech.
,
44
(
3
), pp.
524
531
.10.1016/j.jbiomech.2010.09.010
24.
Zhu
,
Q.
,
Gao
,
X.
,
Li
,
N.
,
Gu
,
W.
,
Eismont
,
F.
, and
Brown
,
M. D.
,
2016
, “
Kinetics of Charged Antibiotic Penetration Into Human Intervertebral Discs: A Numerical Study
,”
J. Biomech.
,
49
(
13
), pp.
3079
3084
.10.1016/j.jbiomech.2016.07.012
25.
Feenstra
,
P. H.
, and
Taylor
,
C. A.
,
2009
, “
Drug Transport in Artery Walls: A Sequential Porohyperelastic-Transport Approach
,”
Comput. Methods Biomech. Biomed. Eng.
,
12
(
3
), pp.
263
276
.10.1080/10255840802459396
26.
Thon
,
M. P.
,
Hemmler
,
A.
,
Glinzer
,
A.
,
Mayr
,
M.
,
Wildgruber
,
M.
,
Zernecke-Madsen
,
A.
, and
Gee
,
M. W.
,
2018
, “
A Multiphysics Approach for Modeling Early Atherosclerosis
,”
Biomech. Model. Mechanobiol.
,
17
(
3
), pp.
617
644
.10.1007/s10237-017-0982-7
27.
Elkin
,
B. S.
,
Shaik
,
M. A.
, and
Morrison
,
B.
,
2010
, “
Fixed Negative Charge and the Donnan Effect: A Description of the Driving Forces Associated With Brain Tissue Swelling and Oedema
,”
Philos. T. Roy. Soc. A
,
368
(
1912
), pp.
585
603
.10.1098/rsta.2009.0223
28.
Chistiakov
,
D. A.
,
Orekhov
,
A. N.
, and
Bobryshev
,
Y. V.
,
2017
, “
Effects of Shear Stress on Endothelial Cells: Go With the Flow
,”
Acta Physiol. (Oxford)
,
219
(
2
), pp.
382
408
.10.1111/apha.12725
29.
Edlich
,
M.
,
Yellowley
,
C. E.
,
Jacobs
,
C. R.
, and
Donahue
,
H. J.
,
2004
, “
Cycle Number and Waveform of Fluid Flow Affect Bovine Articular Chondrocytes
,”
Biorheology
,
41
(
3–4
), pp.
315
–3
22
.https://pubmed.ncbi.nlm.nih.gov/15299264/
30.
Estell
,
E. G.
,
Murphy
,
L. A.
,
Silverstein
,
A. M.
,
Tan
,
A. R.
,
Shah
,
R. P.
,
Ateshian
,
G. A.
, and
Hung
,
C. T.
,
2017
, “
Fibroblast-Like Synoviocyte Mechanosensitivity to Fluid Shear is Modulated by Interleukin-1α
,”
J. Biomech.
,
60
, pp.
91
99
.10.1016/j.jbiomech.2017.06.011
31.
Allen
,
F. D.
,
Hung
,
C. T.
,
Pollack
,
S. R.
, and
Brighton
,
C. T.
,
2000
, “
Serum Modulates the Intracellular Calcium Response of Primary Cultured Bone Cells to Shear Flow
,”
J. Biomech.
,
33
(
12
), pp.
1585
1591
.10.1016/S0021-9290(00)00144-5
32.
Fritton
,
S. P.
, and
Weinbaum
,
S.
,
2009
, “
Fluid and Solute Transport in Bone: Flow-Induced Mechanotransduction
,”
Annu. Rev. Fluid Mech.
,
41
(
1
), pp.
347
374
.10.1146/annurev.fluid.010908.165136
33.
Netz
,
R. R.
,
2020
, “
Mechanisms of Airborne Infection Via Evaporating and Sedimenting Droplets Produced by Speaking
,”
J. Phys. Cem. B
,
124
(
33
), pp.
7093
7101
.10.1021/acs.jpcb.0c05229
34.
Chao
,
C.
,
Wan
,
M.
,
Morawska
,
L.
,
Johnson
,
G.
,
Ristovski
,
Z.
,
Hargreaves
,
M.
,
Mengersen
,
K.
,
Corbett
,
S.
,
Li
,
Y.
,
Xie
,
X.
, and
Katoshevski
,
D.
,
2009
, “
Characterization of Expiration Air Jets and Droplet Size Distributions Immediately at the Mouth Opening
,”
J. Aerosol Sci.
,
40
(
2
), pp.
122
133
.10.1016/j.jaerosci.2008.10.003
35.
Chen
,
C.-C.
, and
Willeke
,
K.
,
1992
, “
Aerosol Penetration Through Surgical Masks
,”
Am. J. Infect. Control
,
20
(
4
), pp.
177
184
.10.1016/S0196-6553(05)80143-9
36.
Lin
,
T.-H.
,
Chen
,
C.-C.
,
Huang
,
S.-H.
,
Kuo
,
C.-W.
,
Lai
,
C.-Y.
, and
Lin
,
W.-Y.
,
2017
, “
Filter Quality of Electret Masks in Filtering 14.6–594 nm Aerosol Particles: Effects of Five Decontamination Methods
,”
PLoS One
,
12
(
10
), p.
e0186217
.10.1371/journal.pone.0186217
37.
Tcharkhtchi
,
A.
,
Abbasnezhad
,
N.
,
Seydani
,
M. Z.
,
Zirak
,
N.
,
Farzaneh
,
S.
, and
Shirinbayan
,
M.
,
2021
, “
An Overview of Filtration Efficiency Through the Masks: Mechanisms of the Aerosols Penetration
,”
Bioactive Mater.
,
6
(
1
), pp.
106
122
.10.1016/j.bioactmat.2020.08.002
38.
Maas
,
S. A.
,
Ellis
,
B. J.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2012
, “
FEBio: Finite Elements for Biomechanics
,”
ASME J. Biomech. Eng.
,
134
(
1
), p.
011005
.10.1115/1.4005694
39.
Ateshian
,
G. A.
,
Nims
,
R. J.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2014
, “
Computational Modeling of Chemical Reactions and Interstitial Growth and Remodeling Involving Charged Solutes and Solid-Bound Molecules
,”
Biomech. Model. Mechanobiol.
,
13
(
5
), pp.
1105
1120
.10.1007/s10237-014-0560-1
40.
Hou
,
J. C.
,
Maas
,
S. A.
,
Weiss
,
J. A.
, and
Ateshian
,
G. A.
,
2018
, “
Finite Element Formulation of Multiphasic Shell Elements for Cell Mechanics Analyses in FEBio
,”
ASME J. Biomech. Eng.
,
140
(
12
), p.
121009
.10.1115/1.4041043
41.
Ateshian
,
G. A.
,
Shim
,
J. J.
,
Maas
,
S. A.
, and
Weiss
,
J. A.
,
2018
, “
Finite Element Framework for Computational Fluid Dynamics in FEBio
,”
ASME J. Biomech. Eng.
,
140
(
2
), p.
021001
.10.1115/1.4038716
42.
Shim
,
J. J.
,
Maas
,
S. A.
,
Weiss
,
J. A.
, and
Ateshian
,
G. A.
,
2019
, “
A Formulation for Fluid–Structure Interactions in FEBio Using Mixture Theory
,”
ASME J. Biomech. Eng.
,
141
(
5
), p.
051010
.10.1115/1.4043031
43.
de Boer
,
R.
,
2000
,
Theory of Porous Media
,
Springer
,
Berlin Heidelberg, Germany
.
44.
Shim
,
J. J.
, and
Ateshian
,
G. A.
,
2021
, “
A Hybrid Biphasic Mixture Formulation for Modeling Dynamics in Porous Deformable Biological Tissues
,”
Arch. Appl. Mech.
, epub.10.1007/s00419-020-01851-8
45.
Shim
,
J. J.
,
Maas
,
S. A.
,
Weiss
,
J. A.
, and
Ateshian
,
G. A.
,
2021
, “
Finite Element Implementation of Biphasic-Fluid Structure Interactions in FEBio
,”
ASME J. Biomech. Eng.
, 143(9), p.
091005
.10.1115/1.4050646
46.
Ateshian
,
G. A.
, and
Ricken
,
T.
,
2010
, “
Multigenerational Interstitial Growth of Biological Tissues
,”
Biomech. Model. Mechanobiol.
,
9
(
6
), pp.
689
702
.10.1007/s10237-010-0205-y
47.
Ateshian
,
G. A.
,
2015
, “
Viscoelasticity Using Reactive Constrained Solid Mixtures
,”
J. Biomech.
,
48
(
6
), pp.
941
947
.10.1016/j.jbiomech.2015.02.019
48.
Zimmerman
,
B. K.
,
Jiang
,
D.
,
Weiss
,
J. A.
,
Timmins
,
L. H.
, and
Ateshian
,
G. A.
,
2021
, “
On the Use of Constrained Reactive Mixtures of Solids to Model Finite Deformation Isothermal Elastoplasticity and Elastoplastic Damage Mechanics
,”
J. Mech. Phys. Solids
,
155
, p.
104534
.10.1016/j.jmps.2021.104534
49.
Bowen
,
R. M.
,
1982
, “
Compressible Porous Media Models by Use of the Theory of Mixtures
,”
Int. J. Eng. Sci.
,
20
(
6
), pp.
697
735
.10.1016/0020-7225(82)90082-9
50.
Bowen
,
R. M.
,
1980
, “
Incompressible Porous Media Models by Use of the Theory of Mixtures
,”
Int. J. Eng. Sci.
,
18
(
9
), pp.
1129
1148
.10.1016/0020-7225(80)90114-7
51.
Holmes
,
M. H.
,
1986
, “
Finite Deformation of Soft Tissue: Analysis of a Mixture Model in Uni-Axial Compression
,”
ASME J. Biomech. Eng.
,
108
(
4
), pp.
372
381
.10.1115/1.3138633
52.
Brinkman
,
H.
,
1949
, “
A Calculation of the Viscous Force Exerted by a Flowing Fluid on a Dense Swarm of Particles
,”
Flow Turbul. Combust.
,
1
(
1
), p.
27
.10.1007/BF02120313
53.
Mauck
,
R. L.
,
Hung
,
C. T.
, and
Ateshian
,
G. A.
,
2003
, “
Modeling of Neutral Solute Transport in a Dynamically Loaded Porous Permeable Gel: Implications for Articular Cartilage Biosynthesis and Tissue Engineering
,”
ASME J. Biomech. Eng.
,
125
(
5
), pp.
602
614
.10.1115/1.1611512
54.
Fick
,
A.
,
1855
, “
V. On Liquid Diffusion
,”
London, Edinburgh, Dublin Philos. Mag. J. Sci.
,
10
(
63
), pp.
30
39
.10.1080/14786445508641925
55.
Kosto
,
K. B.
, and
Deen
,
W. M.
,
2005
, “
Hindered Convection of Macromolecules in Hydrogels
,”
Biophys. J.
,
88
(
1
), pp.
277
86
.10.1529/biophysj.104.050302
56.
Mansour
,
J. M.
, and
Mow
,
V. C.
,
1976
, “
The Permeability of Articular Cartilage Under Compressive Strain and at High Pressures
,”
J. Bone Jt. Surg Am
,
58
(
4
), pp.
509
16
.10.2106/00004623-197658040-00014
57.
Mow
,
V. C.
, and
Mansour
,
J. M.
,
1977
, “
The Nonlinear Interaction Between Cartilage Deformation and Interstitial Fluid Flow
,”
J. Biomech.
,
10
(
1
), pp.
31
39
.10.1016/0021-9290(77)90027-6
58.
Holmes
,
M. H.
,
Lai
,
W. M.
, and
Mow
,
V. C.
,
1985
, “
Singular Perturbation Analysis of the Nonlinear, Flow-Dependent Compressive Stress Relaxation Behavior of Articular Cartilage
,”
ASME J. Biomech. Eng.
,
107
(
3
), pp.
206
218
.10.1115/1.3138545
59.
Holmes
,
M. H.
, and
Mow
,
V. C.
,
1990
, “
The Nonlinear Characteristics of Soft Gels and Hydrated Connective Tissues in Ultrafiltration
,”
J. Biomech.
,
23
(
11
), pp.
1145
1156
.10.1016/0021-9290(90)90007-P
60.
Tinoco
,
I.
,
1995
,
Physical Chemistry: Principles and Applications in Biological Sciences
,
Prentice Hall
,
Englewood Cliffs, NJ
.
61.
Ateshian
,
G. A.
,
Albro
,
M. B.
,
Maas
,
S.
, and
Weiss
,
J. A.
,
2011
, “
Finite Element Implementation of Mechanochemical Phenomena in Neutral Deformable Porous Media Under Finite Deformation
,”
ASME J. Biomech. Eng.
,
133
(
8
), p.
081005
.10.1115/1.4004810
62.
Ogston
,
A. G.
, and
Phelps
,
C. F.
,
1961
, “
The Partition of Solutes Between Buffer Solutions and Solutions Containing Hyaluronic Acid
,”
Biochem. J.
,
78
(
4
), pp.
827
833
.10.1042/bj0780827
63.
Laurent
,
T. C.
, and
Killander
,
J.
,
1964
, “
A Theory of Gel Filtration and Its Exeperimental Verification
,”
J. Chromatogr. A
,
14
, pp.
317
330
.10.1016/S0021-9673(00)86637-6
64.
McNaught
,
A.
,
1997
,
Compendium of Chemical Terminology: IUPAC Recommendations
,
Blackwell Science, Oxford England Malden
,
MA.
65.
Kelly
,
P. D.
,
1964
, “
A Reacting Continuum
,”
Int. J. Eng. Sci.
,
2
(
2
), pp.
129
153
.10.1016/0020-7225(64)90001-1
66.
Eringen
,
A. C.
, and
Ingram
,
J. D.
,
1965
, “
A Continuum Theory of Chemically Reacting Media—I
,”
Int. J. Eng. Sci.
,
3
(
2
), pp.
197
212
.10.1016/0020-7225(65)90044-3
67.
Hou
,
J. S.
,
Holmes
,
M. H.
,
Lai
,
W. M.
, and
Mow
,
V. C.
,
1989
, “
Boundary Conditions at the Cartilage-Synovial Fluid Interface for Joint Lubrication and Theoretical Verifications
,”
ASME J. Biomech. Eng.
,
111
(
1
), pp.
78
87
.10.1115/1.3168343
68.
Overbeek
,
J. T.
,
1956
, “
The Donnan Equilibrium
,”
Prog. Biophys. Biop. Ch.
,
6
, pp.
57
84
.
69.
Katzir-Katchalsky
,
A.
, and
Curran
,
P. F.
,
1965
, “
Nonequilibrium Thermodynamics in Biophysics
,”
Harvard Books in Biophysics
, Vol.
1
,
Harvard University Press
,
Cambridge, UK
.
70.
Beavers
,
G. S.
, and
Joseph
,
D. D.
,
1967
, “
Boundary Conditions at a Naturally Permeable Wall
,”
J. Fluid Mech.
,
30
(
1
), pp.
197
207
.10.1017/S0022112067001375
71.
Ateshian
,
G. A.
,
1991
, “
Biomechanics of Diarthrodial Joints: Applications to the Thumb Carpometacarpal Joint
,” Ph.D. thesis,
Columbia University
,
New York
.