A continuum mechanical framework for modeling tumor growth and treatment in two- and three-phase systems

Abstract

The growth and treatment of tumors is an important problem to society that involves the manifestation of cellular phenomena at length scales on the order of centimeters. Continuum mechanical approaches are being increasingly used to model tumors at the largest length scales of concern. The issue of how to best connect such descriptions to smaller-scale descriptions remains open. We formulate a framework to derive macroscale models of tumor behavior using the thermodynamically constrained averaging theory (TCAT), which provides a firm connection with the microscale and constraints on permissible forms of closure relations. We build on developments in the porous medium mechanics literature to formulate fundamental entropy inequality expressions for a general class of three-phase, compositional models at the macroscale. We use the general framework derived to formulate two classes of models, a two-phase model and a three-phase model. The general TCAT framework derived forms the basis for a wide range of potential models of varying sophistication, which can be derived, approximated, and applied to understand not only tumor growth but also the effectiveness of various treatment modalities.

1 Introduction

Tumor growth and treatment is an area of science of significant interest to society. Ideally, scientists wish to understand fundamental aspects of tumor growth in sufficient detail to enable accurate mathematical models of the behavior at the length scale of interest in humans, which is on the order of centimeters. The problem that arises is a common problem in science and applied mathematics: how to most efficiently and effectively account for important small-scale behavior in larger scale models. To understand the issue more completely, one must consider the scales involved, which we will identify as the molecular scale, the microscale, and the macroscale.

At the molecular scale, one might endeavor, for example, to understand the processes and reactions that lead to the damage and repair of DNA, gene variations important for specific types of cancer, the role of environmental factors, and interactions among contributing factors. Within this context, important fundamental understanding, such as the hallmarks of cancer [24, 25], emerges. Such fundamental, small-scale work is a principal focus of the medical research community and much has been learned over the last few decades. However, a disconnect exists between the molecular scale of such fundamental work and the typical length scale of tumors in humans. As a result, it is not obvious how molecular-scale studies can be used to describe tumor growth [39]. Because of this, tumor growth is often described based upon purely statistical representations of empirical fits to observations [7, 8, 27, 34, 58, 62]. While such fits to data may be good, mechanistic understanding is lacking from such approaches. Put another way, empirical fits are not based on system physics and thus provide an insufficient basis fundamentally to describe factors affecting tumor growth and also to make meaningful, mechanistically based descriptions of how fundamental changes in a system will affect tumor growth.

As an alternative, microscale continuum methods can be used to describe tumor growth. The microscale is a small scale at which continuum mechanical approaches are valid. Even though tumors occur in biological systems, it can be reasoned that common continuum mechanical notions such as conservation principles, mass transfer, reactions, and thermodynamics are applicable and are relatively well understood at the microscale. The challenge that emerges from a microscale modeling approach is the need to represent meaningfully the processes that occur in a wide variety of matter types including healthy tissue, active tumor regions, necrotic tumor regions, extra-cellular matrix regions, and blood vessels. At the microscale, the domains of each of these entities change with time; interfaces form between the entities; and common curves form where three entities meet. Because of this inherent complexity, and the length scales of interest, microscale modeling approaches are not practical or feasible for the mechanistic description of the dynamics of tumors in a living being.

Just as averaging of molecular-scale phenomena is necessary to formulate a microscale continuum model that abstracts away a portion of the mechanistic detail, larger-scale averaging can be performed to derive a macroscale representation from a microscale formulation. At the macroscale, one endeavors to describe the dynamics of the entities (phases, interfaces, common curves) involved in an averaged sense with notions such as volume fractions and other specific entity measures—concepts that do not exist at the microscale. A point in a macroscale model thus represents the averaged conditions embodying all entities around the centroid of a small region. Such models can be formulated in a deterministic sense if and when the averaged conditions are insensitive to small changes in the scale of the averaging region. Useful macroscale models must account for subscale behavior in an approximate, averaged sense; but they must also mechanistically describe tumor dynamics at the scale of applications.

A variety of approaches exist to formulate macroscale models of tumor growth [57]. More broadly considered, a variety of homogenization and averaging methods have been applied to develop macroscale models of porous medium systems [2, 21]. Typically macroscale model formulations are formed and closed phenomenologically directly at the macroscale. Examples for this are the multiparameter models based on mixture theory [26, 37, 50, 51]. While an expedient approach, phenomenological macroscale models do not provide a firm connection with the microscale and cannot be assured to be consistent with the second law of thermodynamics.

Recent continuum mechanics work has resulted in the development of the thermodynamically constrained averaging theory (TCAT) [15, 19, 44]. The TCAT approach formally averages microscale quantities to the macroscale, including not only phases but also interfaces and common curves, incorporates thermodynamics in a scale-consistent manner, and results in entropy inequality expressions that can be used to guide the formulation of models. TCAT also includes evolution equations for the geometry of the phase regions and their boundaries that reduce the closure problem, are based upon mathematical theorems, and are separate from all conservation principles. Recently, notions from integral and differential geometry have been used rigorously to address closure relations needed to describe capillary pressure [42, 47]. Because macroscale TCAT models are firmly connected with microscale antecedents, experimental observations or computational simulations at the microscale can be averaged to the macroscale and used to validate a resultant macroscale model. Considerable microscale resources exist, which can be leveraged to advance macroscale models.

While some aspects of TCAT formulations have been used to model tumor growth [31, 54], a complete and rigorous hierarchy of models formulated and closed using TCAT procedures has not yet been accomplished. Such an advancement is possible by leveraging recent advances in the TCAT approach and applying these to tumor growth and treatment. This advancement provides the development of a hierarchy of models of varying sophistication that can be applied to this important class of problems.

2 Goal and objectives

The overall goal of this work is to develop a framework for modeling tumor growth and treatment at the macroscale. The specific objectives of this work are: (1) to summarize available macroscale conservation, potential, and thermodynamic equations applicable to modeling tumors; (2) to provide a general simplified entropy inequality (SEI), which can be used to guide model closure; (3) to formulate a rigorous two-phase macroscale model for tumor growth; (4) to formulate a rigorous three-phase macroscale model for tumor growth; and (5) to discuss ways in which the model-building framework can be used to formulate models to describe a wide variety of more complex and detailed systems than the ones provided explicitly herein.

3 TCAT approach

The approach to be taken to meet the goal and objectives of this work is to leverage existing TCAT model-building components [19, 52] to formulate models for tumor growth. The advantage of this approach is one of relative simplicity and expediency: using the available formulated components, model building is relatively straightforward; and the substantial amount of effort and manipulations needed to derive an SEI are eliminated. The disadvantage of this approach is that it could be viewed as jumping into the middle of a carefully structured model formulation approach. To circumvent this potential misperception, a brief summary of the TCAT approach is provided to orient readers who are not yet familiar with TCAT. General guidance for the TCAT approach is available in the literature [15, 19, 44,45,46], and specific details of a two-fluid phase compositional model for a porous medium system have also been presented [52].

In general, the TCAT approach is initiated with a general, minimal description of the system to be modeled. This description includes the entities to be modeled, specification of the physical and chemical phenomena occurring within each entity, and the interactions among entities. For the case of concern herein, entities considered include three phases, three interfaces, and a common curve. Entropy, momentum, and energy are resolved at the entity level, and the chemical composition of the mass of each entity is resolved. Classical irreversible thermodynamics is used, and continuum methods are assumed to be valid and deterministic at the macroscale.

All conservation, balance, thermodynamic, and potential equations are formulated at the microscale and then systematically averaged to the macroscale. The macroscale balance of entropy is arranged to solve for the entropy density production rate, which is known to be a nonnegative quantity from the second law of thermodynamics. All macroscale conservation and potential equations are arranged such that the terms in the equation sum to zero. Each collection of terms is multiplied by a Lagrange multiplier and added to a system entropy balance. The Lagrange multipliers are solved for to eliminate material derivatives to the extent possible, yield a dimensionally consistent equation, and connect the processes that produce entropy to the rate of entropy production. Rearrangement of this augmented entropy inequality and reduction to a strict flux–force form, requiring approximations, is a key archival result of the TCAT approach. All of the manipulations leading up to this equation do not need to be repeated for each application that uses the SEI for a given class of models. Furthermore, macroscale conservation and evolution equations are also already available and can be used to formulate models. The chief remaining work when leveraging an extant model hierarchy is thus to use the SEI to formulate model closure relations and combine these equations with conservation and evolution equations to produce a well-posed model. Because of the general approach taken that includes minimal assumptions, a typical SEI supports the formulation of a hierarchy of models of varying sophistication, which results from applying secondary restrictions to the general SEI (e.g. entities of importance and their properties, specific forms of closure relations). This brief overview of the TCAT approach will be detailed in the sections that follow and used to produce example tumor growth models.

We will consider tumor growth models that can be idealized as containing two fluid phases and one solid phase. Compositional effects for mass will be important. An existing TCAT model hierarchy that meets these specifications has been derived [52] and will be relied upon as a foundation for the TCAT modeling approach that follows. Use of this hierarchy will simplify the model formulation process.

4 Macroscale equations

Macroscale equations relied upon in the TCAT approach to form an entropy inequality include conservation equations, a balance of entropy equation, thermodynamic equations, and potential equations. In this section, we summarize only the conservation equations, which are used to construct the target models of concern in this work. Details of the derivation of these equations are available in the literature [19, 52]. The approach followed in this work leverages available results without the burden of reproducing these model components—greatly simplifying the model-building process.

We make use of the conservation equations for mass, momentum, and energy, which can be written in a common form for all species in all entities. The compositional, mass–conservation equation for species ${i}$ in entity ${{\alpha }}$ is

$$\begin{aligned} {\mathcal M}_{**}^{{\overline{\overline{{{i}{{\alpha }}}}}}}&=\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} }\right) }{\mathrm{D}t} + {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}} + \nabla \cdot {\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}}{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}\right) }\nonumber \\&\quad -{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}r^{{i}{{\alpha }}}- {{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}= 0 {\quad \text{ for } {{{i}}\in {{{\mathscr {I}}_{s{}}}}}}, \; {{\alpha }}\in {{\mathscr {I}}}, \end{aligned}$$

(1)

where ${{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}$ is the entity extent measure (volume fraction, specific interfacial area, specific common curve length), ${{\rho ^{{{\alpha }}}}}$ is the mass density, ${{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}$ is the mass fraction, ${\varvec{\mathsf {I}}}$ is the identity tensor, ${{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}}$ is the rate of strain tensor, ${\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}$ is the deviation velocity, ${{r}^{{{{i}{{\alpha }}}}}_{}}$ is the rate of mass production of species ${i}$ resulting from all reactions in entity ${{\alpha }}$, ${\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}$ represents the rate of mass transfer of species ${i}$ from connected entity ${{\kappa }}$ to entity ${{\alpha }}$, ${{\mathscr {I}}_{s{}}}$ is the index set of chemical species, ${{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}$ is the index set of connected entities, and ${{\mathscr {I}}}$ is the index set of all entities in the system. Superscripted entity and species qualifiers denote macroscale quantities.

In general, the set ${{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}$ may contain entities of lower dimensions and higher dimensions than the dimension of the ${{\alpha }}$ entity. For example, the closed set of the solid phase adjoins the wetting fluid phase and the non-wetting fluid phase at the interface that forms between the respective pairs of phases, and the common curve formed at the intersection of the solid phase and the two fluid phases. Thus, ${{{\mathscr {I}}_{\mathrm{c}s}}}=\{{ws},{ns},{wns}\}$, where s in an index specifying the solid phase, and the grouping of indices refer to the respective interface and common curve entities. We will restrict the inter-entity transfer of mass and entropy to entities of at most one dimension higher or lower than a reference entity. For momentum and energy, we will allow a concentrated force along the common curve to act on the solid phase in the most general case, and the inter-entity transfer of internal energy will include interactions between the common curve and the solid phase as well. These restrictions are incorporated into the form of the conservation and balance equations written.

Conservation of momentum may be considered from either a compositional or an overall entity perspective [17]. Taking the latter approach, the momentum equation is written as

$$\begin{aligned} {\varvec{\mathcal P}}_{**}^{{\overline{\overline{{{\alpha }}}}}}&=\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}}\right) }{\mathrm{D}t}+{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}-{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} {\mathbf{g }^{{\overline{{{i}{{\alpha }}}}}}_{}}\nonumber \\&\quad -{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}\left( {\mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}} +{\mathbf{u }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{{i}}}}\right) - {{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}\nonumber \\&\quad -\nabla \cdot \left( {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}}\right) =0 {\quad \text{ for } {{{{\alpha }}}\in {{{\mathscr {I}}}}}}, \end{aligned}$$

(2)

where ${{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}$ is the velocity, ${{{\mathbf{g }^{{\overline{{i}{{\alpha }}}}}_{}}}}$ is the body force per unit mass acceleration vector, $\mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}$ is the velocity of flow in an entity averaged over the boundary of the entity, ${\mathbf{u }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{{i}}}$ is the deviation velocity in an entity averaged over the boundary of the entity, ${{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}}$ is the stress tensor, ${\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}$ represents the transfer of momentum from entity ${{\kappa }}$ to entity ${{\alpha }}$, and singular, or concentrated, forces of a common curve acting on a solid phase are included [19].

Conservation of energy equations is written for an overall entity as

$$\begin{aligned} {\mathcal E}_{**}^{{\overline{\overline{{{\alpha }}}}}}&= \frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}{{{E}^{{\overline{\overline{{{\alpha }}}}}}}}}{\mathrm{D}t} + {{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}\cdot \frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}}\right) }{\mathrm{D}t} + {{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\left( {{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} }\right) \frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} }\right) }{\mathrm{D}t}\nonumber \\&\quad - \frac{{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}\cdot {{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}}{2} {{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} }\right) }{\mathrm{D}t} +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} \frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}}{\mathrm{D}t}{\left( {{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2}}\right) }\nonumber \\&\quad + \left[ {{{{E}^{{\overline{\overline{{{\alpha }}}}}}}}+{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}\displaystyle \frac{{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}\cdot {{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}}{2} +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}\left( {{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2}}\right) }\right] {{\varvec{\mathsf {I}}}} \varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}\nonumber \\&\quad -{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}{\mathbf{g }^{{\overline{{{i}{{\alpha }}}}}}_{}}\cdot \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}+ {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}\right) - {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}h_0^{{\overline{\overline{{{\alpha }}}}}} - {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}h^{{\alpha }}\nonumber \\&\quad -{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}\left[ {{{\overline{E}}_{{i}}^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}} + \frac{\left( {\mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}} + {\mathbf{u }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{{i}}}}\right) \cdot \left( {\mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}} + {\mathbf{u }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{{i}}}}\right) }{2} + {{K}^{\,{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{E{i}}}}\right] \nonumber \\&\quad -{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}\cdot \mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}} -{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{Q_1}}\nonumber \\&\quad -\nabla \cdot \left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}\cdot {{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}+ {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{\mathbf{q }^{{\overline{\overline{{{\alpha }}}}}}}}\right) = 0 {\quad \text{ for } {{{{\alpha }}}\in {{{\mathscr {I}}}}}}, \end{aligned}$$

(3)

where ${{E}^{{\overline{\overline{{{\alpha }}}}}}}$ is the internal energy density, ${{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}$ is the kinetic energy per unit mass resulting from velocity fluctuations, ${\overline{E}}_i^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}$ is the partial mass energy averaged over the boundary of the entity, $K_{Ei}^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}$ is the average deviation kinetic energy averaged over the boundary of the entity, ${\mathbf{q }^{{\overline{\overline{{{\alpha }}}}}}}$ is the non-advective energy flux, $h^{{\alpha }}$ is the energy source density, $h_0^{{\overline{\overline{{{\alpha }}}}}}$ is the energy source density due to body force and velocity fluctuations, and ${\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{Q_1}}$ is the transfer of internal energy from entity ${{\kappa }}$ to entity ${{\alpha }}$ other than by phase change.

Equations (1)–(3) are the basic conservation equations needed to formulate the models of interest in this work. Additional, and available, equations needed for model simplification, closure, and completion will be introduced as needed in the formulation process.

5 Simplified entropy inequality

A key concept from the TCAT approach is the use of an entropy inequality to formulate closure relations that are consistent with the second law of thermodynamics. A strict flux–force form of the entropy inequality is needed to satisfy this purpose; this form is referred to as the SEI. The formulation of an SEI requires skill and substantial mathematical manipulations. However, once derived the SEI can be used without understanding all of the details needed to arrive at the final form. A general SEI is available for the class of model considered in this work [52] and can be written as

$$\begin{aligned}&{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}+{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{p_{}^{{{\alpha }}}}}}{\varvec{\mathsf {I}}}}\right) \varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}+\frac{1}{{{{{\theta }}^{{\overline{\overline{s}}}}}}}\left( {{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{s}}}}} -{{{\epsilon }^{{\overline{\overline{s}}}}}}{\varvec{\mathsf {t}}}^s}\right) \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}}\nonumber \\&\qquad +{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\left[ {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}-{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{\gamma }^{{{\alpha }}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) }\right] \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}}\nonumber \\&\qquad +\frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}\left[ {{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{{\varvec{\mathsf {t}}}^{{\overline{\overline{{wns}}}}}}}+ {{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{{\gamma }^{{wns}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{wns}}}}\right) }\right] \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{wns}}}}}}}\nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}\left[ {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\mathbf{q }^{{\overline{\overline{{{\alpha }}}}}}}+ {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}\mathbf{q }_{\mathbf{g } 0}^{{\overline{\overline{{{\alpha }}}}}} +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}\left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}+ {{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2}}\right) {\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}}\right] \cdot \nabla \left( {\frac{1}{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\right) \nonumber \\&\qquad -{{{{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}\setminus N}\ }\;}} \dfrac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{i\,}{\overline{{{\alpha }}}}}}{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}\cdot \nabla \left[ {{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}+{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} + {{\psi }^{{\overline{{{i}{{\alpha }}}}}}_{}}\right. \nonumber \\&\qquad \left. -\left( {{{\mu }^{{\overline{N{{\alpha }}}}}_{}}+{{K}^{{\overline{\overline{N{{\alpha }}}}}}_E}+\frac{{\mathbf{u }^{{\overline{\overline{N{{\alpha }}}}}}}\cdot {\mathbf{u }^{{\overline{\overline{N{{\alpha }}}}}}}}{2} + {{\psi }^{{\overline{N{{\alpha }}}}}_{}}}\right) \right] \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}+ {{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} + {{{{\psi }^{{\overline{{i}{{\alpha }}}}}_{}}}}}\right) {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}r^{{i}{{\alpha }}}\nonumber \\&\qquad +{{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\langle r_{i\alpha }\psi _{i\alpha }\rangle _{\varOmega _{\alpha },\varOmega }\nonumber \\&\qquad +{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{{i}{{\alpha }}}}}\rightarrow {{i}{{\kappa }}}}{M}} \left[ \frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}} +{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} + {{\psi }^{{\overline{{{i}{{\alpha }}}}}}_{}}}\right) \right. \nonumber \\&\qquad \left. - \frac{1}{{{{\theta }}^{{\overline{\overline{{{\kappa }}}}}}}} \left( {{{\mu }^{{\overline{{i}{{\kappa }}}}}_{}} + {{{{K}^{{\overline{\overline{{i}{{\kappa }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\kappa }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\kappa }}}}}}}}}{2} + {{\psi }^{{\overline{{i}{{\kappa }}}}}_{}}}\right) \right] \nonumber \\&\qquad +{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{c}s}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{i}s}\rightarrow {{i}{{\kappa }}}}{M}} \left[ \frac{1}{{{{{\theta }}^{{\overline{\overline{s}}}}}}}\left( {{{{{\mu }^{{\overline{{i}s}}}_{}}}}+ {{{{K}^{{\overline{\overline{{i}s}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}s}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}s}}}}}}}{2}+ {{{{\psi }^{{\overline{{i}s}}}_{}}}}+\frac{{{\varvec{\mathsf {{\sigma }}}}^{{\overline{\overline{s}}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {C}}}^{{s}}_{}}}{3{{\rho ^{s}}}j^s}}\right) \right. \nonumber \\&\qquad \left. - \frac{1}{{{{\theta }}^{{\overline{\overline{{{\kappa }}}}}}}} \left( {{{\mu }^{{\overline{{i}{{\kappa }}}}}_{}} + {{{{K}^{{\overline{\overline{{i}{{\kappa }}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\kappa }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\kappa }}}}}}}}}{2} + {{\psi }^{{\overline{{i}{{\kappa }}}}}_{}}}\right) \right] \nonumber \\&\qquad +{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{{i}{{\alpha }}}}}\rightarrow {{i}{wns}}}{M}}\Bigg [ \frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}+{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} + {{\psi }^{{\overline{{{i}{{\alpha }}}}}}_{}}}\right) \nonumber \\&\qquad - \frac{1}{{{{\theta }}^{{\overline{\overline{wns}}}}}}\left( {{{\mu }^{{\overline{{i}{wns}}}}_{}}+{{{{K}^{{\overline{\overline{{i}{wns}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{wns}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{wns}}}}}}}}{2} +{{\psi }^{{\overline{{i}{wns}}}}_{}}}\right) \Bigg ] \end{aligned}$$

$$\begin{aligned}&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\Bigg \{{\overset{{{{\alpha }}}\rightarrow {{wn}}}{Q_1}} +{\overset{{{{\alpha }}}\rightarrow {{wn}}}{G_0}} +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {{\overline{E}}_{{i}{{\alpha }}}^{{\overline{{wn}}}} + K_{E{{i}{{\alpha }}}}^{{\overline{\overline{{wn}}}}} +\frac{\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{wn}}}}} \cdot \mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{wn}}}}}}{2} +\psi _{{{i}{{\alpha }}}}^{{\overline{{wn}}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{wn}}}{M}}\nonumber \\&\qquad + \left[ {{\overset{{{{\alpha }}}\rightarrow {{wn}}}{\mathbf{T }_0}}+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{\mathbf{v }^{\,{\overline{\overline{{wn}}}}}_{{{\alpha }}}} -{{\mathbf{v }^{{\overline{s}}}_{}}}}{2}+\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{wn}}}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{wn}}}{M}}}\right] \cdot \left( {{\mathbf{v }^{\,{\overline{\overline{{wn}}}}}_{{{\alpha }}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -\left( {\frac{\mathrm{D}^{{}{{\overline{s}}}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}}{\mathrm{D}t}- \chi _s^{{\overline{\overline{{{\alpha }}s}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t}}\right) p_{{\alpha }}^{wn}\Bigg \}\left( {\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}- \frac{1}{{{{{\theta }}^{{\overline{\overline{{wn}}}}}}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\Bigg \{{\overset{{{{\alpha }}}\rightarrow {{{\alpha }}s}}{Q_1}} + {\overset{{{{\alpha }}}\rightarrow {{{\alpha }}s}}{G_0}} +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {{\overline{E}}_{{i}{{\alpha }}}^{{\overline{{{\alpha }}s}}} + K_{E{{i}{{\alpha }}}}^{{\overline{\overline{{{\alpha }}s}}}} +\frac{\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{{\alpha }}s}}}} \cdot \mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{{\alpha }}s}}}}}{2} +\psi _{{{i}{{\alpha }}}}^{{\overline{{{\alpha }}s}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{{{\alpha }}s}}}{M}} \nonumber \\&\qquad + \left[ {{\overset{{{{\alpha }}}\rightarrow {{{\alpha }}s}}{\mathbf{T }_0}} +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{\mathbf{v }^{\,{\overline{\overline{{{\alpha }}s}}}}_{{{\alpha }}}} -{{\mathbf{v }^{{\overline{s}}}_{}}}}{2}+\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{{\alpha }}s}}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{{{\alpha }}s}}}{M}}}\right] \cdot \left( {{\mathbf{v }^{\,{\overline{\overline{{{\alpha }}s}}}}_{{{\alpha }}}} -{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -\chi _s^{{\overline{\overline{{{\alpha }}s}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t}p_{{\alpha }}^{{{\alpha }}s} \Bigg \}\left( {\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}- \frac{1}{{{{\theta }}^{{\overline{\overline{{{\alpha }}s}}}}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{s}}}^-}\ }\;}} \Bigg \{{\overset{{s}\rightarrow {{{\kappa }}}}{Q_1}} +{\overset{{s}\rightarrow {{{\kappa }}}}{G_0}} + {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {{\overline{E}}_{{i}s}^{{\overline{{{\kappa }}}}} + K_{E{i}s}^{{\overline{\overline{{{\kappa }}}}}} + \frac{\mathbf{u }_{{i}s}^{{\overline{\overline{{{\kappa }}}}}} \cdot \mathbf{u }_{{i}s}^{{\overline{\overline{{{\kappa }}}}}}}{2} + \psi _{{i}s}^{{\overline{{{\kappa }}}}}}\right) {\overset{{{i}s}\rightarrow {{i}{{\kappa }}}}{M}}\nonumber \\&\qquad + \left[ {{\overset{{s}\rightarrow {{{\kappa }}}}{\mathbf{T }_0}}+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{\mathbf{v }^{\,{\overline{\overline{{{\kappa }}}}}}_{s}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}{2} +\mathbf{u }_{{i}s}^{{\overline{\overline{{{\kappa }}}}}}}\right) {\overset{{{i}s}\rightarrow {{i}{{\kappa }}}}{M}}}\right] \cdot \left( {{\mathbf{v }^{\,{\overline{\overline{{{\kappa }}}}}}_{s}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad - \chi _s^{{\overline{\overline{{{\kappa }}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t}\left( {{{\mathbf{n _{s}}\,\cdot \,}}{{\varvec{\mathsf {t}}_{s}}}\cdot {\mathbf{n _{s}}}}\right) _s^{{\kappa }}\Bigg \}\left( {\frac{1}{{{{{\theta }}^{{\overline{\overline{s}}}}}}}- \frac{1}{{{{{\theta }}^{{\overline{\overline{{{\kappa }}}}}}}}}}\right) \nonumber \\&\qquad -\Bigg \{{\overset{{{wn}}\rightarrow {{wns}}}{Q_1}}+ {\overset{{{wn}}\rightarrow {{wns}}}{G_0}}+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {{\overline{E}}_{{i}{wn}}^{{\overline{{wns}}}} + K_{E{i}{wn}}^{{\overline{\overline{{wns}}}}} +\frac{\mathbf{u }_{{i}{wn}}^{{\overline{\overline{{wns}}}}} \cdot \mathbf{u }_{{i}{wn}}^{{\overline{\overline{{wns}}}}}}{2} +\psi _{{i}{wn}}^{{\overline{{wns}}}}}\right) {\overset{{{i}{wn}}\rightarrow {{i}{wns}}}{M}}\nonumber \\&\qquad + \left[ {{\overset{{{wn}}\rightarrow {{wns}}}{\mathbf{T }_0}}+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{wn}}} -{{\mathbf{v }^{{\overline{s}}}_{}}}}{2} + \mathbf{u }_{{i}{wn}}^{{\overline{\overline{{wns}}}}}}\right) {\overset{{{i}{wn}}\rightarrow {{i}{wns}}}{M}}}\right] \cdot \left( {{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{wn}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad + \bigg [\left( {{{{\epsilon }^{{\overline{\overline{{ws}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}}\right) \frac{\mathrm{D}^{{}{{\overline{s}}}}\chi _s^{{\overline{\overline{{ws}}}}}}{\mathrm{D}t} \cos \varphi ^{{\overline{\overline{{ws},{wn}}}}} \nonumber \\&\qquad +\frac{{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}}{{{{\epsilon }^{{\overline{\overline{{ws}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}} \frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \sin \varphi ^{{\overline{\overline{{ws},{wn}}}}}\bigg ]{\gamma }_{wn}^{wns}\Bigg \}\left( {\frac{1}{{{{{\theta }}^{{\overline{\overline{{wn}}}}}}}}- \frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\Bigg \{{\overset{{{{\alpha }}s}\rightarrow {{wns}}}{Q_1}}+ {\overset{{{{\alpha }}s}\rightarrow {{wns}}}{G_0}}+ {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {{\overline{E}}_{{i}{{{\alpha }}s}}^{{\overline{{wns}}}} +K_{E{i}{{{\alpha }}s}}^{{\overline{\overline{{wns}}}}} + \frac{\mathbf{u }_{{i}{{{\alpha }}s}}^{{\overline{\overline{{wns}}}}} \cdot \mathbf{u }_{{i}{{{\alpha }}s}}^{{\overline{\overline{{wns}}}}}}{2} +\psi _{{i}{{{\alpha }}s}}^{{\overline{{wns}}}}}\right) {\overset{{{i}{{\alpha }}s}\rightarrow {{i}{wns}}}{M}}\nonumber \\&\qquad + \left[ {{\overset{{{{\alpha }}s}\rightarrow {{wns}}}{\mathbf{T }_0}}+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{{\alpha }}s}} -{{\mathbf{v }^{{\overline{s}}}_{}}}}{2} + \mathbf{u }_{{i}{{{\alpha }}s}}^{{\overline{\overline{{wns}}}}}}\right) {\overset{{{i}{{\alpha }}s}\rightarrow {{i}{wns}}}{M}}}\right] \cdot \left( {{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{{\alpha }}s}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad + \left( {{{{\epsilon }^{{\overline{\overline{{ws}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}}\right) \frac{\mathrm{D}^{{}{{\overline{s}}}}\chi _s^{{\overline{\overline{{{\alpha }}s}}}}}{\mathrm{D}t} {\gamma }_{{{\alpha }}s}^{wns}\Bigg \}\left( {\frac{1}{{{{\theta }}^{{\overline{\overline{{{\alpha }}s}}}}}} - \frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}}\right) \end{aligned}$$

$$\begin{aligned}&\qquad -{\overset{{s}\rightarrow {{wns}}}{Q_1^*}}\left( {\frac{1}{{{{{\theta }}^{{\overline{\overline{s}}}}}}}-\frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\Bigg \{{{{\eta }^{{\overline{\overline{{{\alpha }}}}}}}}\nabla {{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}-\nabla \left( {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{{p_{}^{{{\alpha }}}}}}\right) \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}\left[ {\nabla \left( {{{{{\mu }^{{\overline{{i}{{\alpha }}}}}_{}}}}+{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} +{{{{\psi }^{{\overline{{i}{{\alpha }}}}}_{}}}}}\right) +{{{\mathbf{g }^{{\overline{{i}{{\alpha }}}}}_{}}}}}\right] \nonumber \\&\qquad -{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}\Bigg [{\overset{{{{\alpha }}}\rightarrow {{{\kappa }}}}{\mathbf{T }_0}}-{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\frac{\left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) }{2}{\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{{\kappa }}}}{M}}\nonumber \\&\qquad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\left( {{\mathbf{v }^{\,{\overline{\overline{{{\kappa }}}}}}_{{{\alpha }}}}+\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{{\kappa }}}}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{{\kappa }}}}{M}}\Bigg ]\Bigg \}\cdot \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}\frac{1}{{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}\Bigg \{ {{{\eta }^{{\overline{\overline{{{\alpha }}}}}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) \cdot \nabla {{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}+ \nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{\gamma }^{{{\alpha }}}}}}\right] \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) \cdot \nabla \left( {{{{{\mu }^{{\overline{{i}{{\alpha }}}}}_{}}}}+{{{{K}^{{\overline{\overline{{i}{{\alpha }}}}}}_E}}}+ \frac{{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}{2} + {{{{\psi }^{{\overline{{i}{{\alpha }}}}}_{}}}}}\right) +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}{{{\mathbf{g }^{{\overline{{i}{{\alpha }}}}}_{}}}}\nonumber \\&\qquad +{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{\alpha }}}^+}\ }\;}}\left[ { {\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}-{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\frac{\left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) }{2}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}+{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\left( {{\mathbf{v }^{\,{\overline{\overline{{{\alpha }}}}}}_{{{\kappa }}}}+\mathbf{u }_{{i}{{\kappa }}}^{{\overline{\overline{{{\alpha }}}}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) {\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}}\right] \nonumber \\&\qquad -\Bigg [{\overset{{{{\alpha }}}\rightarrow {{wns}}}{\mathbf{T }_0}}-{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\frac{\left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) }{2} {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{wns}}}{M}}\nonumber \\&\qquad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\left( {{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{{\alpha }}}}+\mathbf{u }_{{{i}{{\alpha }}}}^{{\overline{\overline{{wns}}}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) {\overset{{{{i}{{\alpha }}}}\rightarrow {{i}{wns}}}{M}}\Bigg ] \Bigg \} \cdot \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -\frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}\Bigg \{{{{\eta }^{{\overline{\overline{{wns}}}}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{wns}}}}\right) \cdot \nabla {{{{\theta }}^{{\overline{\overline{{wns}}}}}}}-\nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{wns}}}}\right) {{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{\gamma }^{{wns}}}}\right] \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{\rho ^{{wns}}}}{{\omega }^{{{i}\,}{\overline{{wns}}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{wns}}}}\right) \cdot \nabla \left( {{{{{\mu }^{{\overline{{i}{wns}}}}_{}}}}+{{{{K}^{{\overline{\overline{{i}{wns}}}}}_E}}}+\frac{{{\mathbf{u }^{{\overline{\overline{{i}{wns}}}}}}}\cdot {{\mathbf{u }^{{\overline{\overline{{i}{wns}}}}}}}}{2}+{{{{\psi }^{{\overline{{i}{wns}}}}_{}}}}}\right) \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{\rho ^{{wns}}}}{{\omega }^{{{i}\,}{\overline{{wns}}}}}{{{\mathbf{g }^{{\overline{{i}{wns}}}}_{}}}}\nonumber \\&\qquad +{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{{wns}}}}^+}\ }\;}}\Bigg [{\overset{{{{\kappa }}}\rightarrow {{wns}}}{\mathbf{T }_0}}-{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\frac{\left( {{{\mathbf{v }^{{\overline{{wns}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) }{2}{\overset{{{i}{{\kappa }}}\rightarrow {{i}{wns}}}{M}}\nonumber \\&\qquad + {{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}\left( {{\mathbf{v }^{\,{\overline{\overline{{wns}}}}}_{{{\kappa }}}}+\mathbf{u }_{{i}{{\kappa }}}^{{\overline{\overline{{wns}}}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) {\overset{{{i}{{\kappa }}}\rightarrow {{i}{wns}}}{M}} \Bigg ]\nonumber \\&\qquad -\langle \mathbf{n} _{\varvec{s}}\cdot \mathbf{t} _{\varvec{s}}^{*}\cdot \mathbf{n} _{\varvec{s}}{} \mathbf{n} _{\varvec{s}}\rangle _{\varOmega _{wns},\varOmega } \Bigg \} \cdot \left( {{{\mathbf{v }^{{\overline{{wns}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad +\frac{1}{{{{{\theta }}^{{\overline{\overline{{wn}}}}}}}}\left[ {\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{w}}}}}}{\mathrm{D}t} - \chi _{s}^{{\overline{\overline{{ws}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} -\frac{{{\gamma }^{{wn}}}{{\hat{k}}}_1^{wn}\left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}-{\epsilon }^{{\overline{\overline{{wn}}}}}_{\mathrm{eq}}}\right) }{p_w^{{wn}}-p_n^{wn}+ {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{wn}}}}{{\omega }^{{{i}\,}{\overline{{wn}}}}}\left( {{{\mathbf{g _{{i}{wn}}}}}\cdot {\mathbf{n _{w}}}}\right) ^{{\overline{{wn}}}}}}\right] \nonumber \\&\qquad \times \left( {p_w^{wn}-p_n^{wn}+ {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{wn}}}}{{\omega }^{{{i}\,}{\overline{{wn}}}}}\left( {{{\mathbf{g _{{i}{wn}}}}}\cdot {\mathbf{n _{w}}}}\right) ^{{\overline{{wn}}}} -{{{\gamma }^{{wn}}}}J_w^{wn}}\right) \end{aligned}$$

$$\begin{aligned}&\qquad +\frac{1}{{{{{\theta }}^{{\overline{\overline{{ws}}}}}}}}\chi _{s}^{{\overline{\overline{{ws}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \left[ {p_w^{ws}+\left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{ws}- {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{ws}}}}{{\omega }^{{{i}\,}{\overline{{ws}}}}}\left( {{{\mathbf{g _{{i}{ws}}}}}\cdot {\mathbf{n _{s}}}}\right) ^{{\overline{{ws}}}} +{{{\gamma }^{{ws}}}}J_s^{ws}}\right] \nonumber \\&\qquad +\frac{1}{{{{{\theta }}^{{\overline{\overline{{ns}}}}}}}}\chi _{s}^{{\overline{\overline{{ns}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \left[ {p_n^{ns}+ \left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{ns}- {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{ns}}}}{{\omega }^{{{i}\,}{\overline{{ns}}}}}\left( {{{\mathbf{g _{{i}{ns}}}}}\cdot {\mathbf{n _{s}}}}\right) ^{{\overline{{ns}}}} +{{{\gamma }^{{ns}}}}J_s^{ns}}\right] \nonumber \\&\qquad -\frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{ws}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}}\right) \frac{\mathrm{D}^{{}{{\overline{s}}}}\chi _s^{{\overline{\overline{{ws}}}}}}{\mathrm{D}t} \bigg [{\gamma }_{wn}^{wns}\cos \varphi ^{{\overline{\overline{{ws},{wn}}}}} +{\gamma }_{ws}^{wns}-{\gamma }_{ns}^{wns}\nonumber \\&\qquad + {{\gamma }^{{wns}}}{\kappa }_{G}^{{\overline{\overline{{wns}}}}} - {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{wns}}}}{{\omega }^{{{i}\,}{\overline{{wns}}}}}\left( {{{\mathbf{g _{{i}{wns}}}}}\cdot \mathbf{n }_{ws}}\right) ^{{\overline{{wns}}}}\bigg ]\nonumber \\&\qquad -\frac{1}{{{{{\theta }}^{{\overline{\overline{{wns}}}}}}}}\left( {\frac{{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}}{{{{\epsilon }^{{\overline{\overline{{ws}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}}}\right) \frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \bigg [{\gamma }_{wn}^{wns}\sin \varphi ^{{\overline{\overline{{ws},{wn}}}}} -{{\gamma }^{{wns}}}{\kappa }_{N}^{{\overline{\overline{{wns}}}}}\nonumber \\&\qquad - \langle \mathbf{n} _{\varvec{s}}\cdot \mathbf{t} _{\varvec{s}}^{*}\cdot \mathbf{n} _{\varvec{s}}\rangle _{\varOmega _{wns},\varOmega _{wns}} + {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\rho ^{{wns}}}}{{\omega }^{{{i}\,}{\overline{{wns}}}}}\left( {{{\mathbf{g _{{i}{wns}}}}}\cdot {\mathbf{n _{s}}}}\right) ^{{\overline{{wns}}}}\bigg ]\nonumber \\&\quad = {{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}{\Lambda }^{{\overline{\overline{{{\alpha }}}}}} \ge 0, \end{aligned}$$

(4)

where all symbols are defined in the notation section and the interactions between the solid phase and the common curve have been explicitly noted.

The general SEI given by Eq. (4) contains the flux–force pairs that can produce entropy in a system. While they may appear overwhelming to the non-specialist, a brief review of the terms by line number can aid understanding of the fluxes considered for which closure relations will be developed. Lines 1–3 are fluxes involving stress tensors and their products with deformation rate tensors. Line 4 is a conductive heat transfer and deviation term flux and their product with a temperature gradient force. Lines 5 and 6 consist of a deviation velocity flux and a force that is a gradient in potential terms, and higher order terms in deviation velocities. Lines 7 and 8 are reaction fluxes resulting from potentials and deviation quantities grouped as forces. Lines 9–14 represent inter-entity mass transfer resulting from differences in potentials and deviation quantities. Lines 18–31 express the inter-entity transfer of energy resulting from forces that are a difference in temperatures. Lines 32–46 are momentum fluxes resulting from differences in entity velocities. Lines 47 and 48 are fluxes in extent measures resulting from a deviation in capillary pressure between the fluid phases from equilibrium. Lines 49 and 50 are changes in porosity resulting from a balance of forces in the direction normal to the solid phase. Lines 51 and 52 represent a change in the wetted fraction of the solid phase resulting from motion of the common curve due to a balance of forces acting tangent to the surface of the solid phase, Lines 53 and 54 express the changes in the porosity resulting from common curve forces acting normal to the solid surface, and Line 55 is the total entropy density production rate of the system. Because each member of the set of fluxes is independent of all other members of the set of fluxes, the fluxes can be considered in turn to derive permissible forms of the closure relations. Examples of how to do so are available in the literature [19, 28, 52]; examples for modeling tumors are detailed in the following sections.

6 Two-phase system

6.1 Description

The purpose of this section is to consider a relatively simple macroscale model to describe tumor growth. This example will enable a straightforward exposition demonstrating the TCAT model building and closure process. Since the focus is on the model formulation process, neither the model solution nor evaluation will be considered herein.

We wish to use the TCAT model components summarized above to formulate a macroscale model consisting of two phases: a solid phase denoted with the index s, and a single fluid phase denoted with the index f, which is a subset of the more general two-fluid phase TCAT model hierarchy that is summarized above. The solid phase consists of the extra-cellular matrix (ECM), live and necrotic tumor cells, glucose, oxygen, water, and a chemotherapeutic drug. These species are important for the solid phase because of the separation of mass transfer and homogeneous phase reactions, and the conceptual representation of the operative processes affecting tumor growth. The fluid phase represents the interstitial fluid, which consists of a dominant water species, glucose, oxygen, and a chemotherapeutic drug. Both phases contain a background species that comprises all other species that are not explicitly considered. Lysis of necrotic cells is also considered.

Based on the above description, the index set of entities is

$$\begin{aligned} {{\mathscr {I}}}=\{f,s,fs\}, \end{aligned}$$

(5)

where fs denotes the fluid–solid interface.

The index set of the nine species considered is denoted

$$\begin{aligned} {{\mathscr {I}}_{s{}}}=\{ l,n,e,g,o,c,w,x,y \}, \end{aligned}$$

(6)

where l denotes a living tumor species, n a necrotic tumor species, e the extra-cellular matrix species, g glucose, o oxygen, c a chemotherapeutic drug, w the water, x a collective background species for the fluid phase, and y a collective background species for the solid phase. The background species represent a set of non-reactive background species that do not undergo mass transfer. The solid phase may contain all species but x, and the fluid phase may contain the g, o, c, w and x species.

Primary restrictions specify the general basis for the TCAT model framework to be used. The primary restrictions for the model are: continuum mechanics represents the system of concern with sufficient fidelity; a clear separation of length scales exists between the microscale and the macroscale; and classical irreversible thermodynamics can be used to describe the system of concern. Secondary restrictions are specified to simplify the general TCAT model hierarchy to the simplest possible form that represents the system of concern. These restrictions have implications for model formulation and closure, including simplification of the general SEI. The secondary restrictions for this application are: (SR1-2P) the system consists of one fluid phase, one solid phase, and an interface; (SR2-2P) the system is isothermal; (SR3-2P) the interface is massless, (SR4-2P) kinetic energy terms are higher order and of negligible importance; (SR5-2P) body force vectors and potentials are identical for all species; (SR6-2P) chemical reactions can be formulated in terms of chemical affinities; (SR7-2P) inertial terms in the momentum equations are insignificant due to the slow dynamics of the systems considered; (SR8-2P) density-weighted, area-averaged velocities and deviation velocities are equal to their volume-averaged counterparts; (SR9-2P) the Lagrangian stress tensor product with the Green’s deformation tensor can be neglected from the driving force difference for mass transfer to the solid phase; and (SR10-2P) the activity coefficient of all species is unity.

6.2 Formulation

The two-phase model described above will be formulated into a closed mathematical model. The steps involved in doing so are detailed in the subsections that follow and include the formulation of a restricted SEI, use of the restricted SEI to formulate a permissible set of closure relations, and a closed model formulation based upon conservation equations and closure approximations.

6.2.1 Restricted SEI

Equation (4) is a general expression that applies to a wide range of systems involving two fluid phases, a solid phase, three interfaces, and a common curve; and a set of flux–force pairs involving dissipative processes involving mass, momentum, and energy. This general expression can be simplified for the example two-phase system described above, and many other applications of concern as well. While several general SEIs have been developed for various systems [16,17,18, 28, 29, 52], the ability to simplify a general SEI for a complex system to describe simpler systems involving fewer entities than the general expression is a hallmark of the TCAT approach. For the application at hand, the original seven entity SEI can be reduced to a three entity SEI, where ${{\mathscr {I}}}=\{f,s,fs\}$. Furthermore, the other secondary restrictions noted above result in substantial simplifications of Eq. (4), while preserving the necessary terms needed to formulate a closed, well-posed model. Applying all secondary restrictions to Eq. (4) yields the restricted SEI given by

$$\begin{aligned}&\frac{1}{{\theta }}\left( {{{\epsilon }^{{\overline{\overline{f}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{f}}}}} +{{\epsilon }^{{\overline{\overline{f}}}}}{{p_{}^{f}}}{\varvec{\mathsf {I}}}}\right) \varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{f}}}}}+\frac{1}{{\theta }}\left( {{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{s}}}}} -{{{\epsilon }^{{\overline{\overline{s}}}}}}{\varvec{\mathsf {t}}}^s}\right) \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}}\nonumber \\&\qquad +\frac{1}{{\theta }}\left[ {{{\epsilon }^{{\overline{\overline{fs}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{fs}}}}} -{{\epsilon }^{{\overline{\overline{fs}}}}}{{\gamma }^{fs}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{fs}}}\right) }\right] \varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{fs}}}}}\nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}\setminus N}\ }\;}} \dfrac{1}{{\theta }}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{i\,}{\overline{{{\alpha }}}}}}{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}\cdot \nabla \left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}- {{\mu }^{{\overline{N{{\alpha }}}}}_{}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}\ }\;}} \frac{1}{{\theta }}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}R^{k{{\alpha }}} A^{k{{\alpha }}} \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{i}f}}\rightarrow {{i}s}}{M}} \frac{1}{{\theta }}\left( {{{\mu }^{{\overline{{i}f}}}_{}} + {{\psi }^{{\overline{f}}}_{}} -{{{{\mu }^{{\overline{{i}s}}}_{}}}}- {{\psi }^{{\overline{s}}}_{}}}\right) \nonumber \\&\qquad -\frac{1}{{\theta }}\Bigg \{-\nabla \left( {{{\epsilon }^{{\overline{\overline{f}}}}} {{p_{}^{f}}}}\right) + {{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\omega }^{{{i}\,}{\overline{f}}}}\nabla {{\mu }^{{\overline{{i}f}}}_{}} +\nabla {{\psi }^{{\overline{f}}}_{}}+{\mathbf{g }^{{\overline{f}}}_{}}}\right) \nonumber \\&\qquad -{\overset{{f}\rightarrow {fs}}{\mathbf{T }_0}} \Bigg \}\cdot \left( {{\mathbf{v }^{{\overline{f}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -\frac{1}{{\theta }}\Bigg \{\nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{fs}}}\right) {{\epsilon }^{{\overline{\overline{fs}}}}}{{\gamma }^{fs}}}\right] +{\overset{{f}\rightarrow {fs}}{\mathbf{T }_0}}+{\overset{{s}\rightarrow {fs}}{\mathbf{T }_0}} \Bigg \} \cdot \left( {{\mathbf{v }^{{\overline{fs}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad +\frac{1}{{\theta }}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \left[ {p_f^{fs} +\left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{fs} +{{\gamma }^{fs}} J_s^{fs}}\right] \nonumber \\&\quad = {{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}{\Lambda }^{{\overline{\overline{{{\alpha }}}}}} \ge 0, \end{aligned}$$

(7)

where ${\theta }$ is the isothermal temperature of the system, ${{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}$ is the index set of all reaction equations in the ${{\alpha }}$ phase, $R^{k{{\alpha }}}$ is a molar rate of reaction, $A^{k{{\alpha }}}$ is a chemical affinity, and SR6-2P has been used to write the reactions in terms of chemical affinities [19]. Each of these flux–force pairs is explained physically and used to formulate a permissible set of closure relations in the section that follows.

6.2.2 Closure relations

We will consider the flux–force pairs in Eq. (7) in order and use these pairs to formulate a permissible set of closure relations, which will be used to formulate a closed, solvable model. References will be made to line numbers in this equation as individual flux–force pairs are considered. Each member of the set of fluxes is independent of all other fluxes, and each member of the set of forces is independent of all other forces. However, a flux may depend upon not only the conjugate force that appears as product in the SEI but also on other members of the set of forces as well, which is referred to as cross-coupled closure. The SEI is an equation that specifies constraints on permissible forms of closure relations. For both of these reasons, closure relations are not unique, and the complexity of the approximate forms can be tailored to the applications. Approaches to derive approximations have been developed and found to have utility for describing a range of physical systems [19]. We will follow these traditional approaches to generate models that can be compared to data. If the models are not consistent with the observations, the closure relation approximations and the model restrictions may be reexamined and modified as needed. Thus, a clear path to modifying macroscale continuum models exists.

The first term in Line 1 of Eq. (7) represents a flux involving the stress tensor for the fluid and the fluid pressure, and the conjugate force is the deformation rate tensor for the fluid. The product of this flux and force must be nonnegative under all conditions and zero at equilibrium. A zero-order approximation is the simplest possible approximation, and it has been found to provide a good description of macroscale porous medium systems [19]. The zero-order approximation is that the flux term is zero under all conditions such that

$$\begin{aligned} {{\varvec{\mathsf {t}}}^{{\overline{\overline{f}}}}}+{{p_{}^{f}}}{\varvec{\mathsf {I}}}=0, \end{aligned}$$

(8)

allowing the macroscale stress tensor to be approximated as

$$\begin{aligned} {{\varvec{\mathsf {t}}}^{{\overline{\overline{f}}}}}=-{{p_{}^{f}}}{\varvec{\mathsf {I}}}. \end{aligned}$$

(9)

This is a reasonable approximation because at the macroscale fluid flow through a porous medium is essentially inviscid, with momentum transfer to the solid phase particles being a dominant process. This dominant exchange process is represented by ${\overset{{f}\rightarrow {fs}}{\mathbf{T }_0}}$, which appears in Line 7 and is discussed below.

Following a similar line of reasoning, the stress tensor for the solid phase is

$$\begin{aligned} {{{\varvec{\mathsf {t}}}^{{\overline{\overline{s}}}}}}=\varvec{\mathsf {t}}^s, \end{aligned}$$

(10)

and the stress tensor for the interface appearing in Line 2 of Eq. (7) is

$$\begin{aligned} {{\varvec{\mathsf {t}}}^{{\overline{\overline{fs}}}}}={{\gamma }^{fs}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{fs}}}\right) . \end{aligned}$$

(11)

Equation (10) suggests that the macroscale stress tensor for the solid phase can be derived based upon intrinsic averaging of the microscale stress tensor, which is determinate based upon the solid behavior at the microscale. Equation (11) relates the interfacial stress tensor to the product of the interfacial tension and the orientation of the interface.

Line 3 in Eq. (7) involves a flux represented by the deviation velocity ${\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}$ and the conjugate force, which is the gradient of the chemical potential. The difference in chemical potential occurs because of the necessary constraint that

$$\begin{aligned} {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}=0 {\quad \text{ for } {{{{\alpha }}}\in {{{\mathscr {I}}}}}}, \end{aligned}$$

(12)

which implies that $N-1$ deviation velocities are independent. Choosing w as the dominant reference species in the f phase and e as the dominant reference species in the s phase yields the linear first-order closure relation for the fluid phase given by

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{f}}}} {\mathbf{u }^{{\overline{\overline{{i}f}}}}}= -{{\hat{\varvec{\mathsf {D}}}}}_{{i}w}^f\cdot \nabla \left( {{{\mu }^{{\overline{{i}f}}}_{}}-{{\mu }^{{\overline{wf}}}_{}}}\right) {\quad \text{ for } {{{i}}\in {\{g,o,c,x \}}}}, \end{aligned}$$

(13)

and the corresponding closure relation for the solid-phase deviation velocity given by

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{s}}}} {\mathbf{u }^{{\overline{\overline{{i}s}}}}}= -{{\hat{\varvec{\mathsf {D}}}}}_{{i}e}^s\cdot \nabla \left( {{{\mu }^{{\overline{{i}s}}}_{}}-{{\mu }^{{\overline{es}}}_{}}}\right) {\quad \text{ for } {{{i}}\in {\{ l,n,g,o,c,w,y\}}}}, \end{aligned}$$

(14)

where ${{\hat{\varvec{\mathsf {D}}}}}_{ij}^{{\alpha }}$ are second-rank symmetric tensors that parameterize the deviation velocities in terms of gradients in chemical potentials.

Line 4 of Eq. (7) expresses the flux–force pair of the rate of reaction and the chemical activity. A permissible closure relation must be generated from this pair and related to the reaction term appearing in Eq. (1), which is the term involving ${{r}^{{{{i}{{\alpha }}}}}_{}}$ that is the rate of mass production of species ${i}$ per time resulting from all reactions in phase ${{\alpha }}$. A linear closure relation is

$$\begin{aligned} R^{k{{\alpha }}}=-{{\hat{K}}}_{k{{\alpha }}}A^{k{{\alpha }}} {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(15)

where ${{\hat{K}}}_{k{{\alpha }}}$ is a nonnegative reaction rate coefficient for the $k^{th}$ reaction in the ${{\alpha }}$ phase. Substitution of this approximation into the SEI yields a nonnegative production rate.

Two tasks remain to formulate closure approximations for ${{r}^{{{{i}{{\alpha }}}}}_{}}$: the closure approximation must be related to the general reaction variable, and the set of chemical reactions occurring in each phase must be formulated. The overall mass production rate results from the set of reactions such that

$$\begin{aligned} {{r}^{{{{i}{{\alpha }}}}}_{}}={{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}\ }\;}} \nu _{{i}k{{\alpha }}}MW_{i}R^{k{{\alpha }}} {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(16)

and the affinity is defined as

$$\begin{aligned} A^{k{{\alpha }}}={{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}\nu _{{i}k{{\alpha }}}MW_{i}{\quad \text{ for } {{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}}}, {{\alpha }}\in {{{\mathscr {I}}_{\mathrm{P}}}}, \end{aligned}$$

(17)

where MW is the molecular weight, and $\nu$ is a molar stoichiometric reaction coefficient. Equations (15)–(17) can be combined to yield the mass production rate resulting from reaction given by

$$\begin{aligned} {{r}^{{{{i}{{\alpha }}}}}_{}}= -{{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}\ }\;}} {{\hat{K}}}_{k{{\alpha }}}\nu _{{i}k{{\alpha }}}MW_{i}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}\nu _{{i}k{{\alpha }}}MW_{i}}\right) {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}. \end{aligned}$$

(18)

If the set of molar biochemical reactions and molecular weights are known for all reactions in each phase, then the species mass production rate is fully specified and closed. When the reactions are nonlinear, resolution issues related to averaging exist, which will not be considered herein [19].

It is assumed that no reactions occur in the fluid phase, which implies that the g, c,  and o species are non-reactive in this phase, although all species may undergo mass transfer to and from the solid phase. Reactions occur in the solid phase that lead to tumor growth, the consumption of glucose, oxygen and the chemotherapeutic drug species, tumor species death due to the chemotherapy drug, conversion of a living tumor to a necrotic tumor resulting from a lack of oxygen, and necrotic tumor lysis. The set of reactions are summarized in Table 1, where $C^{{{i}{{\alpha }}}}$ is the molar concentration of species ${i}$ in phase ${{\alpha }}$, and $MW_{i}$ is the molecular weight of species ${i}$.

Table 1 Reactions for two-phase systems

Line 5 of Eq. (7) is a flux–force pair involving mass transfer between the fluid and solid phases. A conjugate flux–force closure for the mass transfer density rate is

$$\begin{aligned} {\overset{{{i}f}\rightarrow {{i}s}}{M}}= {{\hat{K}}}_M^{{i}fs} \left( {{{\mu }^{{\overline{{i}f}}}_{}}+{{\psi }^{{\overline{f}}}_{}}-{{\mu }^{{\overline{{i}s}}}_{}}-{{\psi }^{{\overline{s}}}_{}}}\right) , \end{aligned}$$

(19)

where ${{\hat{K}}}_M^{{i}fs}$ is a nonnegative mass-transfer rate coefficient.

Lines 6 and 7 of Eq. (7) are a flux–force pair involving momentum of the fluid phase. A conjugate flux–force closure can be formulated for the isotropic case as

$$\begin{aligned}&\nabla \left( {{{\epsilon }^{{\overline{\overline{f}}}}} {{p_{}^{f}}}}\right) -{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{\omega }^{{{i}\,}{\overline{f}}}} \left[ {\nabla \left( {{{\mu }^{{\overline{{i}f}}}_{}} +{{\psi }^{{\overline{f}}}_{}}}\right) +{\mathbf{g }^{{\overline{f}}}_{}}}\right] +{\overset{{f}\rightarrow {fs}}{\mathbf{T }_0}}\nonumber \\&\quad ={{\hat{\varvec{\mathsf {R}}}}}^f\cdot \left( {{\mathbf{v }^{{\overline{f}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) , \end{aligned}$$

(20)

where ${{\hat{\varvec{\mathsf {R}}}}}^f$ is a positive definite resistance tensor for the fluid phase.

Line 8 of Eq. (7) involves the forces acting on the fluid–solid interface, and a conjugate flux–force closure can be written as

$$\begin{aligned} -\nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{fs}}}\right) {{\epsilon }^{{\overline{\overline{fs}}}}}{{\gamma }^{fs}}}\right] -{\overset{{f}\rightarrow {fs}}{\mathbf{T }_0}}-{\overset{{s}\rightarrow {fs}}{\mathbf{T }_0}} = {{\hat{\varvec{\mathsf {R}}}}}^{fs} \cdot \left( {{\mathbf{v }^{{\overline{fs}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) , \end{aligned}$$

(21)

where ${{\hat{\varvec{\mathsf {R}}}}}^{fs}$ is a positive-semidefinite resistance tensor for the interface.

Line 9 of Eq. (7) involves the change in the porosity related to forces acting on the fluid–solid interface. A conjugate flux–force approximation can be written as

$$\begin{aligned} \frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t}= {{\hat{c}}}\left[ {p_f^{fs} +\left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{fs} +{{\gamma }^{fs}} J_s^{fs}}\right] , \end{aligned}$$

(22)

where ${{\hat{c}}}$ is a nonnegative solid compressibility coefficient.

The set of conjugate flux–force closure relations summarized in this subsection can be combined with a set of conservation of mass and momentum equations to yield a closed model. This model is formulated in the section that follows.

6.2.3 Closed model

A closed two-phase model can be formulated by combining conservation equations given in Sect. 4 with closure relations given in Sect. 6.2.2. Because we are assuming isothermal conditions, the model will include conservation of species mass and entity momentum equations, but it will not be necessary to include conservation of energy equations. Other choices of secondary restrictions would result in different macroscale models, but the derivation process would follow a similar path.

The conservation of mass equations for a species in the fluid phase are

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{f}}}}\left( {{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{{\omega }^{{{i}\,}{\overline{f}}}}} }\right) }{\mathrm{D}t} + {{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{{\omega }^{{{i}\,}{\overline{f}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{f}}}}}} + \nabla \cdot \left( {{{\epsilon }^{{\overline{\overline{f}}}}} {\rho ^{f}}{{{\omega }^{{{i}\,}{\overline{f}}}}}{\mathbf{u }^{{\overline{\overline{{i}f}}}}}}\right) \nonumber \\&\quad - {\overset{{{i}s}\rightarrow {{i}f}}{M}} = 0 {\quad \text{ for } {{{i}}\in {\{}}} g,o,c,x \}, \end{aligned}$$

(23)

and for the solid phase

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{s}}}}\left( {{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\rho ^{s}}}{{{\omega }^{{{i}\,}{\overline{s}}}}}}\right) }{\mathrm{D}t} + {{{\epsilon }^{{\overline{\overline{s}}}}}}{{\rho ^{s}}}{{{\omega }^{{{i}\,}{\overline{s}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}} + \nabla \cdot \left( {{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\rho ^{s}}}{{{\omega }^{{{i}\,}{\overline{s}}}}}{\mathbf{u }^{{\overline{\overline{{i}s}}}}}}\right) \nonumber \\&\quad -{{{\epsilon }^{{\overline{\overline{s}}}}}}r^{{i}s} - {\overset{{{i}f}\rightarrow {{i}s}}{M}} = 0 {\quad \text{ for } {{{i}}\in {\{}}} l,n,g,o,c,w,y \}, \end{aligned}$$

(24)

where the dominant species for each phase is eliminated using the constraint equation

$$\begin{aligned} {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}=1 {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(25)

where ${{{\mathscr {I}}_{\mathrm{P}}}}$ is the index set of phases and some species are omitted from each phase by definition. Thus, four species conservation of mass equations exist for the fluid phase, and seven species conservation of mass equations exist for the solid phase. Combining these equations with constraints given by Eq. (25) fully specifies the composition of both phases.

Equations (23) and (24) contain terms involving deviation velocities, reactions, and interphase mass transfer. The previously derived closure relations can be used to explicitly specify these equations. Substituting Equations (14) and (19) into Eq. (23) yields the species conservation of mass equation for the fluid phase

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{f}}}}\left( {{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{{\omega }^{{{i}\,}{\overline{f}}}}} }\right) }{\mathrm{D}t} + {{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{{\omega }^{{{i}\,}{\overline{f}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{f}}}}}} - \nabla \cdot \left[ {{{\epsilon }^{{\overline{\overline{f}}}}} {\rho ^{f}} {{\hat{\varvec{\mathsf {D}}}}}_{{i}w}^f\left( {{{\mu }^{{\overline{if}}}_{}}-{{\mu }^{{\overline{wf}}}_{}}}\right) }\right] \nonumber \\&\quad + {{\hat{K}}}_M^{{i}fs} \left( {{{\mu }^{{\overline{{i}f}}}_{}}+{{\psi }^{{\overline{f}}}_{}}-{{\mu }^{{\overline{{i}s}}}_{}}-{{\psi }^{{\overline{s}}}_{}}}\right) = 0 {\quad \text{ for } {{{i}}\in {\{}}} g,o,c,x \}, \end{aligned}$$

(26)

and for the solid phase

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{s}}}}\left( {{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}}{{{\omega }^{{{i}\,}{\overline{s}}}}} }\right) }{\mathrm{D}t} + {{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}}{{{\omega }^{{{i}\,}{\overline{s}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}} - \nabla \cdot \left[ {{{\epsilon }^{{\overline{\overline{f}}}}} {\rho ^{f}} {{\hat{\varvec{\mathsf {D}}}}}_{{i}e}^s\left( {{{\mu }^{{\overline{is}}}_{}}-{{\mu }^{{\overline{es}}}_{}}}\right) }\right] \nonumber \\&\quad +{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{i}s}}}\ }\;}} {{\hat{K}}}_{ks}\nu _{{i}k s}MW_{i}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\mu }^{{\overline{{i}s}}}_{}}\nu _{{i}k s}MW_{i}}\right) \nonumber \\&\quad - {{\hat{K}}}_M^{{i}fs} \left( {{{\mu }^{{\overline{{i}f}}}_{}}+{{\psi }^{{\overline{f}}}_{}}-{{\mu }^{{\overline{{i}s}}}_{}}-{{\psi }^{{\overline{s}}}_{}}}\right) = 0{\quad \text{ for } {{{i}}\in {\{}}} l,n,g,o,c,w,y \}, \end{aligned}$$

(27)

where only the reactions involving species ${i}$ need to be considered for the ${i}$ species transport equation, which is specified by the summation over index set ${{\mathscr {I}}_{\mathrm{rxn}{i}s}}$ that is the set of all reactions in the s entity that involve species ${i}$.

Chemical potentials can be written in terms of mass fractions to minimize the closure problem with the conservation of mass equations. The macroscale chemical potential may be written as

$$\begin{aligned} {{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}= \mu _0^{{\overline{{{i}{{\alpha }}}}}}({{{p_{}^{{{\alpha }}}}}},{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}) +\frac{R_g{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}}}{MW_{i}}\ln \left( {{{x}^{{\overline{\overline{{{i}{{\alpha }}}}}}}}\hat{\gamma }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}\right) , \end{aligned}$$

(28)

where $\mu _0^{{\overline{{{i}{{\alpha }}}}}}({{{p_{}^{{{\alpha }}}}}},{{{{\theta }}^{{\overline{\overline{{{\alpha }}}}}}}})$ is a reference state chemical potential, $R_g$ is the ideal gas constant, ${{x}^{{\overline{\overline{{{i}{{\alpha }}}}}}}}$ is a mole fraction, and $\hat{\gamma }^{{\overline{\overline{{{i}{{\alpha }}}}}}}$ is an activity coefficient, which we will assume is equal to 1 for all species and all compositions encountered. The mole fraction can be written as

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}=\frac{MW_{i}}{MW_{{\alpha }}}{{x}^{{\overline{\overline{{{i}{{\alpha }}}}}}}}, \end{aligned}$$

(29)

where the molecular weight of the ${{\alpha }}$ entity is defined as

$$\begin{aligned} MW_{{\alpha }}={{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}MW_{i}{{x}^{{\overline{\overline{{{i}{{\alpha }}}}}}}}={{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}\left( {\frac{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}}{MW_{i}}}\right) ^{-1}. \end{aligned}$$

(30)

Equation (1) can also be summed over all species, which yields an overall conservation of mass equation for an entity. This equation can be further expanded using an approximation for the macroscale velocity, which follows from the conservation of momentum equations that follow. These standard manipulations are detailed in the literature [19].

Recall the conservation of momentum equation, which can be written as

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}}\right) }{\mathrm{D}t} +{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}-{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{\mathbf{g }^{{\overline{{{\alpha }}}}}_{}}\nonumber \\&\quad -{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}\left( {\mathbf{v }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}} +{\mathbf{u }^{{\overline{\overline{{{\alpha }},{{\kappa }}}}}}_{{i}}}}\right) - {{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}\nonumber \\&\quad -\nabla \cdot \left( {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}}\right) =0 {\quad \text{ for } {{{{\alpha }}}\in {{{\mathscr {I}}}}}}. \end{aligned}$$

(31)

This equation may be combined with the closure relations and manipulated to deduce a pair of closed momentum equations for the phases. The interface momentum equation is trivial because of the secondary restriction specifying a massless interface. Equations (9), (19), and (20) can be substituted into Eq. (31) to yield the resultant fluid phase momentum equation

$$\begin{aligned}&{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{{\omega }^{{{i}\,}{\overline{f}}}} \nabla {{\mu }^{{\overline{{i}f}}}_{}} -{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{\mathbf{g }^{{\overline{f}}}_{}} +{{\hat{\varvec{\mathsf {R}}}}}^f\cdot \left( {{\mathbf{v }^{{\overline{f}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\quad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{\hat{K}}}_M^{{i}fs} \left( {{{\mu }^{{\overline{{i}f}}}_{}}+{{\psi }^{{\overline{f}}}_{}}-{{\mu }^{{\overline{{i}s}}}_{}}-{{\psi }^{{\overline{s}}}_{}}}\right) \left( {{\mathbf{v }^{{\overline{f}}}_{}}+{\mathbf{u }^{{\overline{\overline{{i}f}}}}}}\right) =0, \end{aligned}$$

(32)

or approximated using the Gibbs–Duhem equation as [19]

$$\begin{aligned}&{{\epsilon }^{{\overline{\overline{f}}}}}\nabla {{p_{}^{f}}} -{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{\mathbf{g }^{{\overline{f}}}_{}} +{{\hat{\varvec{\mathsf {R}}}}}^f\cdot \left( {{\mathbf{v }^{{\overline{f}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\quad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{\hat{K}}}_M^{{i}fs} \left( {{{\mu }^{{\overline{{i}f}}}_{}}+{{\psi }^{{\overline{f}}}_{}}-{{\mu }^{{\overline{{i}s}}}_{}}-{{\psi }^{{\overline{s}}}_{}}}\right) \left( {{\mathbf{v }^{{\overline{f}}}_{}}+{\mathbf{u }^{{\overline{\overline{{i}f}}}}}}\right) =0, \end{aligned}$$

(33)

where the inertial terms have been dropped due to the slow dynamics of the system and the area averaged velocities multiplying the mass exchange approximations have been replaced with the volume averages, as specified in secondary restrictions SR7-2P and SR8-2P.

Equation (31) can be summed over all entities, eliminating the momentum exchange terms, and yielding after substitution of the closure relations for the macroscale stress tensors

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{f}}}}\left( {{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}} {\mathbf{v }^{{\overline{f}}}_{}}}\right) }{\mathrm{D}t}+\frac{\mathrm{D}^{{}{{\overline{s}}}}\left( {{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}} {\mathbf{v }^{{\overline{s}}}_{}}}\right) }{\mathrm{D}t} +{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}} {\mathbf{v }^{{\overline{f}}}_{}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{f}}}}} +{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}} {\mathbf{v }^{{\overline{s}}}_{}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}\nonumber \\&\quad -{{\epsilon }^{{\overline{\overline{f}}}}}{\rho ^{f}}{\mathbf{g }^{{\overline{f}}}_{}} -{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}}{\mathbf{g }^{{\overline{s}}}_{}} -\nabla \cdot \left[ {-{{\epsilon }^{{\overline{\overline{f}}}}}{{p_{}^{f}}}{\varvec{\mathsf {I}}}+ {{{\epsilon }^{{\overline{\overline{s}}}}}}\varvec{\mathsf {t}}^s +{{\gamma }^{fs}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{fs}}}\right) }\right] =0. \end{aligned}$$

(34)

Further manipulations using this equation are possible, and expressions for the stress tensor have been developed using TCAT [20]. The end objective of such a momentum equation is to determine the velocity of the solid phase. Conditions may exist in which this velocity can be approximated without consideration of the details of solid mechanics. The existence of mass exchange complicates the usual solid mechanics mapping between Eulerian and Lagrangian coordinate systems, and open research questions remain, which will not be considered in this work.

7 Three-phase system

7.1 Description

For the two-phase system considered above, the tumor species were considered part of the solid phase. However, it is generally accepted that cells and tissues may be modeled as fluids [14]. Such approaches may be used to model adhesion of cells among themselves, to the ECM, and to substrates [1, 30, 56], and to model chemotaxis and haptotaxis to represent active cellular movement, such as observed in angiogenesis, invasion, and branching [9, 31].

Hence, the purpose of this section is to consider a more complicated three-phase model involving tumor growth. This will require defining the entities, species in the entities, reactions, and other aspects of the model. An alternative form of the SEI will be relied upon to guide model closure. The basic model-building steps will parallel those steps followed to formulate the two-phase model presented in the previous section.

The three phases consist of a solid phase, s, a wetting fluid denoted w, and a non-wetting fluid phase denoted n. The wetting fluid preferentially wets the solid phase compared to the non-wetting fluid phase, although a contact angle can exist between the fluid–fluid interface and the solid surface; measured through the wetting fluid phase this contact angle will be $<90$ degrees and 0 degrees for a strongly water wet solid phase. The index set of all entities is

$$\begin{aligned} {{\mathscr {I}}}=\{w,n,s,wn, ws, ns, wns\}, \end{aligned}$$

(35)

where the grouping of two symbols denotes the interface that forms between the two respective phases and wns is an index for the common curve formed where three phases meet.

The index set of the ten species considered is identical to the set previously considered for the two-phase tumor model with the exception that an additional collective background species z is added to the set

$$\begin{aligned} {{\mathscr {I}}_{s{}}}=\{ l,n,e,g,o,c,w,x,y,z \}. \end{aligned}$$

(36)

Table 2 Composition of three-phase systems

Table 2 summarizes the entities, the index used for each entity, and species composition of each entity in the system; the composition of interfaces and common curves will not be resolved in this restricted model, so the composition of these entities is neglected. The primary restrictions for the three-phase model are identical to the primary restrictions for the previously presented two-phase model. The secondary restrictions are similar to the secondary restrictions for the two-phase model, although not identical: (SR1-3P) the system is isothermal; (SR2-3P) the interfaces and common curves are massless, (SR3-3P) kinetic energy terms are higher order and of negligible importance; (SR4-3P) body force vectors and potentials are identical for all species; (SR5-3P) chemical reactions can be formulated in terms of chemical affinities; (SR6-3P) inertial terms in the momentum equations are insignificant due to the slow dynamics of the systems considered; (SR7-3P) density-weighted, area-averaged velocities and deviation velocities are equal to their volume-averaged counterparts; (SR8-3P) the Lagrangian stress tensor product with the Green’s deformation tensor can be neglected from the driving force difference for mass transfer to the solid phase; (SR9-3P) the activity coefficient of all species is unity; (SR10-3P) common curve effects are considered higher order and neglected; (SR11-3P) the contact angle is considered constant as a result of the relatively slow dynamics of tumor growth; and (SR12-3P) the curvature of the solid phase is not dependent upon the fluid wetting the solid surface.

7.2 Formulation

The three-phase model described above is formulated into a closed model following a similar approach that was taken for the simpler two-phase model formulated above. The sections that follow detail the restricted SEI, the closure relations, and the closed model.

7.2.1 Restricted SEI

Applying the secondary restrictions specified in Sect. 7.1 to the general SEI given by Eq. (4) results in the restricted SEI for a compositional three-phase system given by

$$\begin{aligned}&{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}\frac{1}{{\theta }}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}+{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{p_{}^{{{\alpha }}}}}}{\varvec{\mathsf {I}}}}\right) \varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}+\frac{1}{{\theta }}\left( {{{{\epsilon }^{{\overline{\overline{s}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{s}}}}} -{{{\epsilon }^{{\overline{\overline{s}}}}}}{\varvec{\mathsf {t}}}^s}\right) \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}}}\nonumber \\&\qquad +{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}\frac{1}{{\theta }}\left[ {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}-{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{\gamma }^{{{\alpha }}}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) }\right] \varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}}\nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}\setminus N}\ }\;}} \dfrac{1}{{\theta }}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{\omega }^{{i\,}{\overline{{{\alpha }}}}}}{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}\cdot \nabla \left( {{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}- {{\mu }^{{\overline{N{{\alpha }}}}}_{}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}\ }\;}} \frac{1}{{\theta }}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}R^{k{{\alpha }}} A^{k{{\alpha }}}\nonumber \\&+{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{i}w}}\rightarrow {{i}n}}{M}} \frac{1}{{\theta }}\left( {{{\mu }^{{\overline{{i}w}}}_{}} +{{\psi }^{{\overline{w}}}_{}} -{{\mu }^{{\overline{{i}n}}}_{}} - {{\psi }^{{\overline{n}}}_{}}}\right) \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{i}w}}\rightarrow {{i}s}}{M}} \frac{1}{{\theta }}\left( {{{\mu }^{{\overline{{i}w}}}_{}} + {{\psi }^{{\overline{w}}}_{}} -{{\mu }^{{\overline{{i}s}}}_{}} - {{\psi }^{{\overline{s}}}_{}}}\right) \nonumber \\&\qquad +{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{\overset{{{{i}n}}\rightarrow {{i}s}}{M}} \frac{1}{{\theta }}\left( {{{\mu }^{{\overline{{i}n}}}_{}} + {{\psi }^{{\overline{n}}}_{}} -{{\mu }^{{\overline{{i}s}}}_{}} - {{\psi }^{{\overline{s}}}_{}}}\right) \nonumber \\&\qquad - \frac{1}{{\theta }}\Bigg \{-\nabla \left( {{{\epsilon }^{{\overline{\overline{w}}}}} {{p_{}^{w}}}}\right) + {{{\epsilon }^{{\overline{\overline{w}}}}}}{{\rho ^{w}}}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\omega }^{{{i}\,}{\overline{w}}}}\nabla {{\mu }^{{\overline{{i}w}}}_{}}+\nabla {{{\psi }^{{\overline{w}}}_{}}}+{{\mathbf{g }^{{\overline{w}}}_{}}}}\right) \nonumber \\&\qquad -{\overset{{w}\rightarrow {{wn}}}{\mathbf{T }_0}}-{\overset{{w}\rightarrow {{ws}}}{\mathbf{T }_0}}\Bigg \}\cdot \left( {{{\mathbf{v }^{{\overline{w}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad - \frac{1}{{\theta }}\Bigg \{-\nabla \left( {{{\epsilon }^{{\overline{\overline{n}}}}} {{p_{}^{n}}}}\right) + {{{\epsilon }^{{\overline{\overline{n}}}}}}{{\rho ^{n}}}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\omega }^{{{i}\,}{\overline{n}}}}\nabla {{\mu }^{{\overline{{i}n}}}_{}}+\nabla {{{\psi }^{{\overline{n}}}_{}}}+{{\mathbf{g }^{{\overline{n}}}_{}}}}\right) \nonumber \\&\qquad -{\overset{{n}\rightarrow {{wn}}}{\mathbf{T }_0}}-{\overset{{n}\rightarrow {{ns}}}{\mathbf{T }_0}}\Bigg \}\cdot \left( {{{\mathbf{v }^{{\overline{n}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad -{{\sum \limits _{{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}\frac{1}{{\theta }}\Bigg \{\nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{{\gamma }^{{{\alpha }}}}}}\right] +{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{\alpha }}}^+}\ }\;}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}\Bigg \} \cdot \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\qquad +\frac{1}{{\theta }}\left[ {\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{w}}}}}}{\mathrm{D}t} - \chi _{s}^{{\overline{\overline{{ws}}}}} \frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} -\frac{{{\gamma }^{{wn}}}{{\hat{k}}}_1^{wn}\left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}-{\epsilon }^{{\overline{\overline{{wn}}}}}_{\mathrm{eq}}}\right) }{p_w^{{wn}}-p_n^{wn}}}\right] \nonumber \\&\qquad \times \left( {p_w^{wn}-p_n^{wn}-{{{\gamma }^{{wn}}}}J_w^{wn}}\right) \nonumber \\&\qquad +\frac{1}{{\theta }}\chi _{s}^{{\overline{\overline{{ws}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \left[ {p_w^{ws}+\left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{ws}+{{{\gamma }^{{ws}}}}J_s^{ws}}\right] \nonumber \\&\qquad +\frac{1}{{\theta }}\chi _{s}^{{\overline{\overline{{ns}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} \left[ {p_n^{ns}+ \left( {{{\mathbf{n _{s}}\,\cdot \,}}\varvec{\mathsf {t}}_s \cdot {\mathbf{n _{s}}}}\right) _s^{ns}+{{{\gamma }^{{ns}}}}J_s^{ns}}\right] \nonumber \\&\quad = {{\sum \limits _{{{{\alpha }}}\in {{{\mathscr {I}}}}\ }\;}}{\Lambda }^{{\overline{\overline{{{\alpha }}}}}} \ge 0. \end{aligned}$$

(37)

Equation (37) can be used to generate closure relations that are consistent with the second law of thermodynamics. While the form of this equation is much simpler than the general case given by Eq. (4) because of the secondary restrictions applied, it is significantly more complicated than the two-phase SEI given by Eq. (7) as a result of the additional entities that must be considered.

7.2.2 Closure relations

The procedure suggested for deriving closure relations is similar to that used for the two-phase case, although additional relations are needed for the three-phase case. We consider the first 16 lines of Eq. (37) to generate a candidate set of closure relations, which will in turn be used to formulate a closed model in the following section.

A zero-order closure for the fluid phases yields from Line 1 of Eq. (37)

$$\begin{aligned} {{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}}=-{{{p_{}^{{{\alpha }}}}}}{\varvec{\mathsf {I}}}{\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}}}, \end{aligned}$$

(38)

and for the solid phase

$$\begin{aligned} {{{\varvec{\mathsf {t}}}^{{\overline{\overline{s}}}}}}=\varvec{\mathsf {t}}^s. \end{aligned}$$

(39)

A zero-order approximation based on Line 2 yields an approximation for the stress tensor for the interfaces

$$\begin{aligned} {{{\varvec{\mathsf {t}}}^{{\overline{\overline{{{\alpha }}}}}}}}={{\gamma }^{{{\alpha }}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}}}. \end{aligned}$$

(40)

A first-order conjugate flux–force approximation based upon Line 3 in Eq. (37) yields the following approximations for the deviation velocities in each of the phases

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{w}}}} {\mathbf{u }^{{\overline{\overline{{i}w}}}}}&= -{{\hat{\varvec{\mathsf {D}}}}}_{{i}w}^w\cdot \nabla \left( {{{\mu }^{{\overline{{i}w}}}_{}}-{{\mu }^{{\overline{ww}}}_{}}}\right) {\quad \text{ for } {{{i}}\in {\{g,o,c,z \}}}}, \end{aligned}$$

(41)

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{n}}}} {\mathbf{u }^{{\overline{\overline{{i}n}}}}}&= -{{\hat{\varvec{\mathsf {D}}}}}_{{i}l}^n\cdot \nabla \left( {{{\mu }^{{\overline{{i}n}}}_{}}-{{\mu }^{{\overline{ln}}}_{}}}\right) {\quad \text{ for } {{{i}}\in {\{n,g,o,c,w,x \}}}}, \end{aligned}$$

(42)

and

$$\begin{aligned} {{\omega }^{{{i}\,}{\overline{s}}}} {\mathbf{u }^{{\overline{\overline{{i}s}}}}}= -{{\hat{\varvec{\mathsf {D}}}}}_{{i}e}^s\cdot \nabla \left( {{{\mu }^{{\overline{{i}s}}}_{}}-{{\mu }^{{\overline{es}}}_{}}}\right) {\quad \text{ for } {{{i}}\in {\{g,o,w,y \}}}}, \end{aligned}$$

(43)

where water is taken as the reference species in the wetting phase, the live tumor species is the reference species in the non-wetting phase, and the extra-cellular matrix species is the reference species for the solid phase.

Table 3 Reactions for three-phase systems

Line 4 in Eq. (37) can be used in a conjugate flux–force form to approximate thermodynamically consistent reaction rates. The manner in which this is done is analogous to the approach previously detailed for the simpler two-phase model. An example set of reactions are detailed in Table 3. The reactions are related to tumor growth, destruction, necrotic species formation, and lysis. All reactions occur in the non-wetting fluid phase. We emphasize that this is merely an example set of reactions and alternative reaction sets are permissible. The general expression for reactions in phases is

$$\begin{aligned} {{r}^{{{{i}{{\alpha }}}}}_{}}=- {{\sum \limits _{{k}\in {{{\mathscr {I}}_{\mathrm{rxn}{{\alpha }}}}}\ }\;}} {{\hat{K}}}_{k{{\alpha }}}\nu _{{i}k{{\alpha }}}MW_{i}\left( {{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{\mu }^{{\overline{{{i}{{\alpha }}}}}}_{}}\nu _{{i}k{{\alpha }}}MW_{i}}\right) {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}. \end{aligned}$$

(44)

A first-order conjugate form closure relation for the rate of mass transfer can be deduced from Lines 5–7 in Eq. (37) and is of the form

$$\begin{aligned} {\overset{{{i}{\alpha }}\rightarrow {{i}{\beta }}}{M}}= {{\hat{K}}}_M^{{i}{\alpha \beta }} \left( {{{\mu }^{{\overline{{i}{\alpha }}}}_{}}+{{\psi }^{{\overline{{\alpha }}}}_{}}-{{\mu }^{{\overline{{i}{\beta }}}}_{}}-{{\psi }^{{\overline{{\beta }}}}_{}}}\right) {\quad \text{ for } {{{\alpha },{\beta }}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(45)

where the indices denote that mass transfer can occur between any of three binary combinations of phases.

Lines 8–11 of Eq. (37) can be used to generate cross-coupled approximations involving the interphase transfer of momentum for the fluid phases of the form

$$\begin{aligned}&\nabla \left( {{{{\epsilon }^{{\overline{\overline{w}}}}}}{{{p_{}^{w}}}}}\right) -{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{w}}}}}}{{\rho ^{w}}}{{\omega }^{{{i}\,}{\overline{w}}}} \left[ {\nabla \left( {{{\mu }^{{\overline{{i}w}}}_{}} +{{{\psi }^{{\overline{w}}}_{}}}}\right) +{{\mathbf{g }^{{\overline{w}}}_{}}}}\right] +{\overset{{w}\rightarrow {{wn}}}{\mathbf{T }_0}}+{\overset{{w}\rightarrow {{ws}}}{\mathbf{T }_0}}\nonumber \\&\quad ={{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}{{\hat{\varvec{\mathsf {R}}}}}^w_{{\kappa }}\cdot \left( {{\mathbf{v }^{{\overline{{{\kappa }}}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) , \end{aligned}$$

(46)

and

$$\begin{aligned}&\nabla \left( {{{{\epsilon }^{{\overline{\overline{n}}}}}}{{{p_{}^{n}}}}}\right) -{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{n}}}}}}{{\rho ^{n}}}{{\omega }^{{{i}\,}{\overline{n}}}} \left[ {\nabla \left( {{{\mu }^{{\overline{{i}n}}}_{}} +{{{\psi }^{{\overline{n}}}_{}}}}\right) +{{\mathbf{g }^{{\overline{n}}}_{}}}}\right] +{\overset{{n}\rightarrow {{wn}}}{\mathbf{T }_0}}+{\overset{{n}\rightarrow {{ns}}}{\mathbf{T }_0}}\nonumber \\&\quad ={{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}{{\hat{\varvec{\mathsf {R}}}}}^n_{{\kappa }}\cdot \left( {{\mathbf{v }^{{\overline{{{\kappa }}}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) , \end{aligned}$$

(47)

where ${{\hat{\varvec{\mathsf {R}}}}}$ is a positive semi-definite resistance tensor.

A first-order approximation based upon Line 12 in Eq. (37) yields

$$\begin{aligned}&\nabla \cdot \left[ {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\alpha }}}}}\right) {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\gamma }^{{{\alpha }}}}}\right] +{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{\alpha }}}^+}\ }\;}}{\overset{{{{\kappa }}}\rightarrow {{{\alpha }}}}{\mathbf{T }_0}}\nonumber \\&\quad =-{{\hat{\varvec{\mathsf {R}}}}}_{{\alpha }}^{{\alpha }}\cdot \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) -{{\sum \limits _{{{{\kappa }}}\in {{{\mathscr {I}}}_{\mathrm{c}{{\alpha }}}^+}\ }\;}}{{\hat{R}}}_{{\kappa }}^{{\alpha }}\cdot \left( {{\mathbf{v }^{{\overline{{{\kappa }}}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}}}. \end{aligned}$$

(48)

Lines 13 and 14 from Eq. (37) can be used to develop a first-order approximation for the relaxation of capillary pressure to its equilibrium state, which is the case when the term in Line 14 vanishes. The resulting approximation is

$$\begin{aligned} \frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{w}}}}}}{\mathrm{D}t} - \chi _{s}^{{\overline{\overline{{ws}}}}} \frac{\mathrm{D}^{{}{{\overline{s}}}}\epsilon }{\mathrm{D}t} -\frac{{{\gamma }^{{wn}}}{{\hat{k}}}_1^{wn}\left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}- {\epsilon }^{{\overline{\overline{{wn}}}}}_{\mathrm{eq}}}\right) }{p_w^{{wn}}-p_n^{wn}} = {{\hat{c}}}^{wn}\left( {p_w^{wn}-p_n^{wn}-{{{\gamma }^{{wn}}}}J_w^{wn}}\right) , \end{aligned}$$

(49)

where ${{\hat{c}}}^{wn}$ is a positive capillary pressure relaxation coefficient.

As a result of SR12-3P, the curvature of the solid phase simplifies to

$$\begin{aligned} J_s^{ss}=J_s^{ws}=J_s^{ns}, \end{aligned}$$

(50)

which provide an approximation for the change in porosity as a function of normal forces on the solid surface of the form

$$\begin{aligned} \frac{\mathrm{D}^{{}{{\overline{s}}}}{\epsilon }}{\mathrm{D}t}= {{\hat{c}}}^{ss} \left[ {\chi _{s}^{{\overline{\overline{{ws}}}}}p_w^{wn}+\chi _{s}^{{\overline{\overline{{ns}}}}}p_n^{ns}+\left( {{{\mathbf{n _{s}}\,\cdot \,}}{{\varvec{\mathsf {t}}_{s}}}\cdot {\mathbf{n _{s}}}}\right) ^{ss} +\left( {\chi _{s}^{{\overline{\overline{{ws}}}}}{{\gamma }^{{ws}}}+\chi _{s}^{{\overline{\overline{{ns}}}}}{{\gamma }^{{ns}}}}\right) J_s^{ss}}\right] , \end{aligned}$$

(51)

where ${{\hat{c}}}^{ss}$ is a positive compressibility coefficient.

In addition to the closure relations that can be deduced from the SEI, evolution equations based upon averaging theorems [19, 22] can also be used to close the three-phase model. The use of such evolution equations to produce closed models is an important aspect of TCAT models. For the three-phase model considered here, an evolution equation for the fluid–fluid interfacial area can be written as

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{{wn}}}}}}}{\mathrm{D}t} + \nabla \cdot \left[ {{{\epsilon }^{{\overline{\overline{{wn}}}}}}\left( {\mathbf{w }^{wn}-{\varvec{\mathsf {G}}^{{wn}}}\cdot {\mathbf{v }^{{\overline{s}}}_{}}}\right) }\right] + {{\epsilon }^{{\overline{\overline{{wn}}}}}}{\varvec{\mathsf {G}}^{{wn}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}} - J_w^{wn}\left( {\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{w}}}}}}{\mathrm{D}t} +\chi ^{{\overline{\overline{{ws}}}}}_s \frac{\mathrm{D}^{{}{{\overline{s}}}}{{{\epsilon }^{{\overline{\overline{s}}}}}}}{\mathrm{D}t}}\right) \nonumber \\&\qquad - {{\hat{k}}}^{wn}\left( {{\epsilon }^{{\overline{\overline{{wn}}}}}_\mathrm{{eq}}-{{\epsilon }^{{\overline{\overline{{wn}}}}}}}\right) - \cos {\varphi ^{{\overline{\overline{ws,wn}}}}}({{\epsilon }^{{\overline{\overline{{ws}}}}}}+{{\epsilon }^{{\overline{\overline{{ns}}}}}}) \frac{\mathrm{D}^{{}{{\overline{s}}}}\chi _s^{{\overline{\overline{{ws}}}}}}{\mathrm{D}t} \nonumber \\&\qquad + \sin {\varphi ^{{\overline{\overline{ws,wn}}}}}\frac{{{\epsilon }^{{\overline{\overline{{wns}}}}}}}{{{\epsilon }^{{\overline{\overline{{ws}}}}}}+{{\epsilon }^{{\overline{\overline{{ns}}}}}}}\frac{\mathrm{D}^{{}{{\overline{s}}}}{{{\epsilon }^{{\overline{\overline{s}}}}}}}{\mathrm{D}t}=0, \end{aligned}$$

(52)

where neglecting common curve contributions this equation simplifies to

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{{wn}}}}}}}{\mathrm{D}t} + \nabla \cdot \left[ {{{\epsilon }^{{\overline{\overline{{wn}}}}}}\left( {\mathbf{w }^{wn}-{\varvec{\mathsf {G}}^{{wn}}}\cdot {\mathbf{v }^{{\overline{s}}}_{}}}\right) }\right] + {{\epsilon }^{{\overline{\overline{{wn}}}}}}{\varvec{\mathsf {G}}^{{wn}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}} - J_w^{wn}\left( {\frac{\mathrm{D}^{{}{{\overline{s}}}}{{\epsilon }^{{\overline{\overline{w}}}}}}{\mathrm{D}t} +\chi ^{{\overline{\overline{{ws}}}}}_s \frac{\mathrm{D}^{{}{{\overline{s}}}}{{{\epsilon }^{{\overline{\overline{s}}}}}}}{\mathrm{D}t}}\right) \nonumber \\&\qquad - {{\hat{k}}}^{wn}\left( {{\epsilon }^{{\overline{\overline{{wn}}}}}_\mathrm{{eq}}-{{\epsilon }^{{\overline{\overline{{wn}}}}}}}\right) - \cos {\varphi ^{{\overline{\overline{ws,wn}}}}}({{\epsilon }^{{\overline{\overline{{ws}}}}}}+{{\epsilon }^{{\overline{\overline{{ns}}}}}}) \frac{\mathrm{D}^{{}{{\overline{s}}}}\chi _s^{{\overline{\overline{{ws}}}}}}{\mathrm{D}t} =0. \end{aligned}$$

(53)

${\varvec{\mathsf {G}}^{{wn}}}$ is a diagonal tensor of trace 1 for the isotropic case, which is a reasonable approximation. For other cases, an evolution equation for this quantity is needed, which can be derived at the microscale, averaged to the macroscale, and approximated in a convenient form.

The interfacial velocity vector is given by $\mathbf{w }^{wn}$, which is the macroscale velocity of the fluid–fluid interface, which moves in a direction normal to the interface. An exact equation for this kinematic quantity is not available. An approximation can be derived based upon the averaging theorems of the form

$$\begin{aligned} {{{\epsilon }^{{\overline{\overline{{wn}}}}}}}^2\mathbf{w }^{wn}= -\frac{\partial {{{{\epsilon }^{{\overline{\overline{w}}}}}}}}{\partial {t}}\left( {\nabla {{{\epsilon }^{{\overline{\overline{w}}}}}}+\langle \mathbf{n} _{w}\rangle _{\varOmega _{ws},\varOmega }}\right) , \end{aligned}$$

(54)

where an approximation is needed for the average of the normal vector, which is trivial for the common isotropic case.

As surfaces deform, their curvatures change. An evolution equation relating the change in the mean and the Gaussian curvature can be written as

$$\begin{aligned}&{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}\frac{\partial {K_n^{wn}}}{\partial {t}} + K_n^{wn}\frac{\partial {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}}}{\partial {t}}+\nabla \cdot \left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}K_n^{wn}\mathbf{w }^{wn}}\right) \nonumber \\&\quad +\frac{J_n^{wn}}{2}\nabla \nabla \varvec{\mathsf {:}}\left[ {\left( {{\varvec{\mathsf {I}}}-{{\varvec{\mathsf {G}}}^{{{wn}}}_{}}}\right) \frac{\partial {{{{\epsilon }^{{\overline{\overline{n}}}}}}}}{\partial {t}}}\right] \nonumber \\&\quad +\frac{J_n^{wn}}{2}\nabla \cdot \left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}J_n^{wn}\mathbf{w }^{wn}}\right) =0. \end{aligned}$$

(55)

A constraint equation based upon Galilean invariance for the mean curvature is

$$\begin{aligned}&-\nabla \cdot \left[ {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}J_n^{wn}\left( {{\varvec{\mathsf {I}}}-2{\varvec{\mathsf {G}}^{{wn}}}}\right) }\right] +2K_n^{wn}\nabla {{{\epsilon }^{{\overline{\overline{n}}}}}}\nonumber \\&\quad -\nabla \nabla \varvec{\mathsf {:}}\left( {\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{wn}}}}\right) \nabla {{{\epsilon }^{{\overline{\overline{n}}}}}}}\right) =0, \end{aligned}$$

(56)

and the Gaussian curvatures can be related to the Euler characteristic for smooth closed boundaries using the Gauss Bonnet theorem giving

$$\begin{aligned} 4\pi {\chi }^{{\overline{\overline{n}}}}= {{{\epsilon }^{{\overline{\overline{{wn}}}}}}}K_n^{wn}+ {{\epsilon }^{{\overline{\overline{{ns}}}}}}K_n^{ns}. \end{aligned}$$

(57)

Additional closure relations for the three-phase model are available through the specification of state equations. Another example of a state equation is a formulation that can be deduced from integral geometry that relates volume fractions, interfacial areas, and curvatures of the form [47]

$$\begin{aligned} {{{\epsilon }^{{\overline{\overline{n}}}}}}= F\left( {{{{\epsilon }^{{\overline{\overline{{wn}}}}}}}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}, {{{\epsilon }^{{\overline{\overline{{wn}}}}}}}J_w^{wn}+{{{\epsilon }^{{\overline{\overline{{ns}}}}}}}J_s^{ns},{\chi }^{{\overline{\overline{n}}}}}\right) , \end{aligned}$$

(58)

where F is a smooth differentiable function that can be deduced from state data. Equation (58) has been shown to provide an accurate representation of capillary pressure that is based on theory, hysteretic free—unlike traditional capillary pressure equations, and applicable for not only equilibrium conditions but also dynamic states as well [23]. Equations of state can be formulated to relate densities to pressures and compositions as well.

7.2.3 Closed model

The purpose of this section is to assemble the pieces needed for a complete, closed three-phase model of tumor growth, which can be used as a basis to derive numerical approximations. The three-phase model is more complicated than the previously formulated two-phase model. This results from the existence of additional entities and the attendant closure problem. Overall conservation equations for mass and momentum are needed for each phase, and compositional equations including mass transfer and reactions are needed for the species in a phase. The massless interface and common curve assumptions do simplify the formulation for lower dimensional entities. We note that other TCAT formulation approaches are possible, and we are merely providing an example for completeness. This model follows from the conservation equations and closure relations based upon the SEI that were previously presented. The quantities we wish to resolve include the composition of each phase, the velocity and mass density of each phase, fluid pressures and solid stresses, entity extents, and interfacial curvatures needed to approximate the state of the system. The model that follows is based upon these considerations.

Summing Eq. (1) over all species yields a set of conservation of mass equations for the phases of the form

$$\begin{aligned} \frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}}\right) }{\mathrm{D}t} + {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}} - {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}=0 {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(59)

which can be written using Eq. (45) as

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}}\right) }{\mathrm{D}t} + {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}} - {{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{{\hat{K}}}_M^{{i}{{\kappa }}{{\alpha }}} \left( {{{\mu }^{{\overline{{i}{{\kappa }}}}}_{}}+{{\psi }^{{\overline{{{\kappa }}}}}_{}}-{{\mu }^{{\overline{{i}{{\alpha }}}}}_{}}-{{\psi }^{{\overline{{{\alpha }}}}}_{}}}\right) =0 \nonumber \\&\qquad {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}, \end{aligned}$$

(60)

where the chemical potentials can be written in terms of mole fractions and mass fractions as previously formulated in Eqs. (28) and (30).

Species composition for each phase follows from

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{{{\alpha }}}}}}\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} }\right) }{\mathrm{D}t} + {{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}} {\varvec{\mathsf {I}}}\varvec{\mathsf {:}}{{{\varvec{\mathsf {d}}}^{{\overline{\overline{{{\alpha }}}}}}}} + \nabla \cdot {\left( {{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}{\rho ^{{{\alpha }}}}}{{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}}{{\mathbf{u }^{{\overline{\overline{{i}{{\alpha }}}}}}}}}\right) }\nonumber \\&\quad -{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}r^{{i}{{\alpha }}}- {{{{\sum \limits _{{{{\kappa }}}\in {{{{{\mathscr {I}}_{\mathrm{c}{{\alpha }}}}}}}\ }\;}}}}{\overset{{{i}{{\kappa }}}\rightarrow {{{i}{{\alpha }}}}}{M}}= 0 {\quad \text{ for } {{{i}}\in {{{\mathscr {I}}_{s{}}}}}}, \; {{\alpha }}\in {{{\mathscr {I}}_{\mathrm{P}}}}, \end{aligned}$$

(61)

where the deviation velocities follow from substitution of Eqs. (41)–(43), reaction terms from Eq. (44), and mass transfer from Eq. (45). Note that a complete set of mass conservation equations must be formulated and solved subject to the constraint

$$\begin{aligned} {{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}=1 {\quad \text{ for } {{{{\alpha }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}}}. \end{aligned}$$

(62)

Conservation of mass equations are not required for interfaces or the common curve because they have been specified as massless.

As a result of secondary restrictions SR6-3P and SR7-3P, Eq. (2) can be simplified, Eqs. (38), (45), and (46) substituted in and the resulting conservation of momentum equations for the fluid phases written as

$$\begin{aligned}&{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\omega }^{{{i}\,}{\overline{{{\alpha }}}}}}\nabla {{\mu }^{{\overline{{i}{{\alpha }}}}}_{}} -{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{g }^{{\overline{{{\alpha }}}}}_{}}}+{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}{{\hat{\varvec{\mathsf {R}}}}}_{{\kappa }}^{{\alpha }}\cdot \left( {{\mathbf{v }^{{\overline{{{\kappa }}}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\quad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\hat{K}}}_M^{{i}{{\alpha }}{{\kappa }}} \left( {{{\mu }^{{\overline{{i}{{\alpha }}}}}_{}}+{{{\psi }^{{\overline{{{\alpha }}}}}_{}}}-{{\mu }^{{\overline{{i}{{\kappa }}}}}_{}}-{{\psi }^{{\overline{{{\kappa }}}}}_{}}}\right) \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}+{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}}\right) =0 , \end{aligned}$$

(63)

which can be approximated using the Gibbs–Duhem equation as [19]

$$\begin{aligned}&{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}\nabla {{{p_{}^{{{\alpha }}}}}}-{{{\epsilon }^{{\overline{\overline{{{\alpha }}}}}}}}{{\rho ^{{{\alpha }}}}}{{\mathbf{g }^{{\overline{{{\alpha }}}}}_{}}}+{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{f}}}}}\ }\;}}{{\hat{\varvec{\mathsf {R}}}}}_{{\kappa }}^{{\alpha }}\cdot \left( {{\mathbf{v }^{{\overline{{{\kappa }}}}}_{}}-{{\mathbf{v }^{{\overline{s}}}_{}}}}\right) \nonumber \\&\quad +{{{{\sum \limits _{{{i}}\in {{{\mathscr {I}}_{s{}}}}\ }\;}}}}{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{P}}}}}\ }\;}}{{\hat{K}}}_M^{{i}{{\alpha }}{{\kappa }}} \left( {{{\mu }^{{\overline{{i}{{\alpha }}}}}_{}}+{{{\psi }^{{\overline{{{\alpha }}}}}_{}}}-{{\mu }^{{\overline{{i}{{\kappa }}}}}_{}}-{{\psi }^{{\overline{{{\kappa }}}}}_{}}}\right) \left( {{{\mathbf{v }^{{\overline{{{\alpha }}}}}_{}}}+{\mathbf{u }^{{\overline{\overline{{{i}{{\alpha }}}}}}}}}\right) =0 , \end{aligned}$$

(64)

where viscous coupling of momentum transport between the fluid phases follows naturally from the cross-coupled closure relations [36].

A momentum equation is also required for the solid phase, which can alternatively be expressed as the total change in momentum for the system, which can be formulated summing Eq. (2) over all entities using the closure approximations given by Eqs. (38) and (40) yielding

$$\begin{aligned}&\frac{\mathrm{D}^{{}{{\overline{w}}}}\left( {{{\epsilon }^{{\overline{\overline{w}}}}}{\rho ^{w}} {\mathbf{v }^{{\overline{w}}}_{}}}\right) }{\mathrm{D}t} +\frac{\mathrm{D}^{{}{{\overline{n}}}}\left( {{{\epsilon }^{{\overline{\overline{n}}}}}{\rho ^{n}} {\mathbf{v }^{{\overline{n}}}_{}}}\right) }{\mathrm{D}t} +\frac{\mathrm{D}^{{}{{\overline{s}}}}\left( {{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}} {\mathbf{v }^{{\overline{s}}}_{}}}\right) }{\mathrm{D}t} +{{\epsilon }^{{\overline{\overline{w}}}}}{\rho ^{w}} {\mathbf{v }^{{\overline{w}}}_{}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{w}}}}} +{{\epsilon }^{{\overline{\overline{n}}}}}{\rho ^{n}} {\mathbf{v }^{{\overline{n}}}_{}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{n}}}}}\nonumber \\&\quad +{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}} {\mathbf{v }^{{\overline{s}}}_{}}{{\varvec{\mathsf {I}}}}\varvec{\mathsf {:}}{{\varvec{\mathsf {d}}}^{{\overline{\overline{s}}}}} -{{\epsilon }^{{\overline{\overline{w}}}}}{\rho ^{w}}{\mathbf{g }^{{\overline{w}}}_{}} -{{\epsilon }^{{\overline{\overline{n}}}}}{\rho ^{n}}{\mathbf{g }^{{\overline{n}}}_{}} -{{\epsilon }^{{\overline{\overline{s}}}}}{\rho ^{s}}{\mathbf{g }^{{\overline{s}}}_{}}\nonumber \\&\quad -\nabla \cdot \left[ {-{{\epsilon }^{{\overline{\overline{w}}}}}{{p_{}^{w}}}{\varvec{\mathsf {I}}}-{{\epsilon }^{{\overline{\overline{n}}}}}{{p_{}^{n}}}{\varvec{\mathsf {I}}}+ {{{\epsilon }^{{\overline{\overline{s}}}}}}\varvec{\mathsf {t}}^s +{{\sum \limits _{{{{\kappa }}}\in {{{{\mathscr {I}}_{\mathrm{I}}}}}\ }\;}}{{\gamma }^{{{\kappa }}}}\left( {{\varvec{\mathsf {I}}}-{\varvec{\mathsf {G}}^{{{\kappa }}}}}\right) }\right] =0 . \end{aligned}$$

(65)

In the event that a zero velocity at the macroscale is not a reasonable assumption for the solid phase, expressions available in the literature [20] can be used to produce a solvable equation. As previously noted, additional work on the solid mechanics of these biomechanical systems is warranted, and this is an active area of research [6, 11].

The three-phase model is complete when augmented with evolution equations for the volume fraction given by Eq. (49) and the interfacial area given by Eq. (53), and equations of state given by Eq. (58) and state equations for the dependency of phase densities as a function of pressures and compositions. Approximations would also be needed for the specific interfacial area ${{\epsilon }^{{\overline{\overline{ns}}}}}$ and the mean curvature of the solid phase $J_s^{ns}$.

8 Discussion

Tumor occurrence and growth exhibits classical problems of scale. A complete understanding requires molecular, genetic, and cellular perspectives, while the systems of concern are at scales of the human body and thus several orders of magnitude above the length scales of fundamental concern. The perspective advanced herein is a practical approach focused on advancing a mechanistic mathematical description that is able to describe the system of concern. The TCAT approach provides a means to formulate such models, which are founded upon conservation principles, thermodynamics, and a set of mathematical theorems. Even though molecular and cellular approaches are essential to advance fundamental understanding of tumor formation and growth, the laws of continuum mechanics must still apply provided the necessary primary restrictions have been met. The approach taken is to develop macroscale models in which tumors are described in an averaged sense with admissible changes in both space and time. With this general approach, a wide variety of specific models is possible.

The primary goal of this work was to illustrate how the TCAT approach can be used to formulate macroscale models of varying sophistication and complexity. Two specific examples were provided: a two-phase formulation, and three-phase formulation; the former included a single fluid phase and the latter included two fluid phases and both included a solid phase. Both of these examples included interfaces, and the full three-phase case can include a common curve, which was included in the SEI but dropped for simplicity in the final example formulation. Both approaches include several species, reactions for tumor formation, tumor destruction, necrotic tissue formation, and necrotic lysis, and mass transfer. We emphasize that the examples provided are intended to be a starting point for advancement in understanding of the how the TCAT approach can be applied, evaluated, and validated. It is expected that details of the methods, especially the species and reaction sets, will evolve as understanding of the essential components of an optimally useful model also evolve. As understanding, and available computing power, expands, so too will the complexity and fidelity of the underlying models. However, the approaches illustrated will provide a solid foundation for continuum mechanical modeling of tumor growth.

Some previous tumor growth models have included elements of TCAT but not to the extent detailed in this work [32, 53,54,55, 60]. They have been successful in modeling tumor growth but not entirely satisfactory. Especially the missing interfaces required questionable approximations when linking capillary pressures with saturations. It is hoped that the presented approach will allow a step forward in this direction. The next steps beyond the numerical implementation of the detailed equations will be the inclusion of angiogenesis in the model and finally the addition of drug delivery. The number of parameters involved will require a sensitivity analysis to identify the leading order effects and to motivate the most important experimental work. Uncertainty quantification is strongly linked to this problem. It is hoped that mechanistic models of the sort presented herein will help to elucidate mechanisms resulting in cancer growth and treatment, which may be considered independent of the genetic origins. In fact, recent studies have revealed extensive variations in genetic signatures, gene expression, and post-translational modifications among different tumors and within tumors, creating complex tumor heterogeneities. Heterogeneity is frequently attributed to genetics, but it is related to non-genetic influences too [5]. Clearly, this directly impacts therapeutic treatments and outcomes [40, 43, 59]. Thus, better understanding of heterogeneities will help to improve treatments of metastatic diseases. This brings us to the final aim of our modeling effort which is the evaluation of the efficiency of cancer drugs. This needs an efficient model of tumor growth linked to a bio-distribution model of the drug under scrutiny.

A motivation for three-fluid-phase models would be a need to represent tumor systems based upon three distinct phases with immiscible fluid-like properties that form separate regions of matter with distinct interfacial tensions. While such models may be needed, much can be done with the general framework for one-fluid and two-fluid models advanced in this work. In the event that model evaluation and validation efforts of the sorts of models advanced in this work result in a conclusion that three-fluid-phase continuum models are needed to represent tumor growth with adequate fidelity, TCAT can be used to derive such models. While the conservation and balance equations are in place for such advancements, additional thermodynamic, evolution equation, and state equation work would be required. A new general SEI would also need to be derived. Of course, such a theoretical framework would be applicable to others systems as well.

Models with some similarities to those detailed herein are being used for successful medical applications (e.g. [6, 38]), and extensions to clinical applications of drug delivery and evaluation of drug efficiency are also possible; initial works in these directions can be found in [4, 10].

Because the intent of this work was to provide a foundational example of how TCAT can be used to formulate macroscale models to describe tumor growth, it has not included the numerical methods aspects of approximating these models. Furthermore, the worth of a model is dependent upon its ability to represent physical systems of concern with sufficient fidelity to be of use to those assessing such systems. Such an assessment requires not only an approximation of the formulation but also parameter estimation and comparison to observed systems. It is hoped that this work will enable efforts to approximate, apply, evaluate, and validate a range of candidate TCAT models in pursuit of realistic and worthwhile representations of systems important to society. Much work remains to be done to fulfill this vision.

9 Conclusions

Several conclusions follow from this work:

• 1. Macroscale continuum mechanical approaches provide a means to describe tumor formation, growth, and various treatment modalities at a length scale relevant to human systems of concern.

• 2. TCAT provides a framework for formulating such macroscale models in a manner that is consistent across length scales with respect to established conservation and thermodynamic principles.

• 3. Two-phase and three-phase TCAT models were formulated as examples of the formulation process and reduced to a closed form.

• 4. The work advances prior efforts in providing a more complete set of entities and including recent advances in differential geometry evolution equations and integral geometry state equations, which provide accurate closure approximation for two-fluid porous medium models.

• 5. The foundational nature of this work is intended to support future work to improve upon the example formulations provided in this work in concert with numerical approximation of the formulated models, parameter estimation, evaluation, and validation.

• 6. It is anticipated that advances in fundamental molecular and cellular level understanding of tumor formation, growth, and treatment processes will inform and help improve continuum scale models of focus in this work, which may enable more routine simulation at scales of concern and help advance more effective treatment approaches.