P.K. Sahoo et al

*Corresponding author name,

# Energy Conditions in Non-minimally Coupled $f(R,T)$ Gravity

P.K. Sahoo  Sanjay Mandal  Simran Arora  Department of Mathematics, Birla
Institute of Technology and Science-Pilani, Hyderabad Campus,
Hyderabad-500078, India<EMAIL_ADDRESS>
(22 Sept 2020; 10 Oct 2020; 19 Oct 2020)

###### Abstract

In today’s scenario, going beyond Einstein’s theory of gravity leads us to
some more complete and modified gravity theories. One of them is the $f(R,T)$
gravity in which $R$ is the Ricci scalar, and $T$ is the trace of the energy-
momentum tensor. Using a well-motivated linear $f(R,T)$ gravity model with a
single parameter, we studied the strong energy condition (SEC), the weak
energy condition (WEC), the null energy condition (NEC), and the dominant
energy condition (DEC) under the simplest non-minimal matter geometry coupling
with a perfect fluid distribution. The model parameter is constrained by
energy conditions and a single parameter proposed equation of state (EoS),
resulting in the compatibility of the $f(R,T)$ models with the accelerated
expansion of the universe. It is seen that the EoS parameter illustrate the
quintessence phase in a dominated accelerated phase, pinpoint to the
cosmological constant yields as a prediction the phantom era. Also, the
present values of the cosmological constant and the acceleration of the
universe are used to check the viability of our linear $f(R,T)$ model of
gravity. It is observed that the positive behavior of DEC and WEC indicates
the validation of the model. In contrast, SEC is violating the condition
resulting in the accelerated expansion of the universe.

###### keywords:

cosmology: theory, cosmological parameters, large-scale structure, cosmology:
observations

††articletype: Article Type

##  1 Introduction

Verification of the universe late-time acceleration has led to comprehensive
studies into its interpretation. Observational data on type Ia supernovae
revealed an acceleration of the current universe (Riess et al., 1998, 1999;
Perlmutter et al., 1999). This phenomenon can either be overruled by
introducing a new kind of universal fluid known as dark energy or proposing
modification and variations of the Einstein theory of gravity. This is done by
modifying the left-hand side of the Einstein equation, i.e., the spacetime’s
geometry. Modified theories of gravity are the geometrical generalizations of
General Relativity (GR) in which the Einstein-Hilbert action is modified by
replacing the Ricci scalar R with a more general function along with matter-
geometry coupling. Some widely known and used modified theories are $f(R)$
gravity (Nojiri Odinstov, 2003), $f(T)$ gravity (Femaro Fiorini, 2007), $f(G)$
gravity (Nojiri ., 2008), $f(R,T)$ gravity (Harko et al, 2011), $f(Q,T)$
gravity (Yixin et al, 2019) etc. In the last few years, $f(R,T)$ gravity
gained much attention where $R$ is the Ricci scalar, and $T$ is the trace of
the energy-momentum tensor.

A lot of research on various aspects of $f(R,T)$ gravity such as
thermodynamics (Sharif Zubair, 2019), scalar pertubations (Alvarenga et al.,
2013), wormhole solutions (Moraes Sahoo, 20171) has been discussed in
literature. Apart from these effective experiments, any modified theory of
gravity should follow an important role in the energy conditions. These
conditions define the attractive existence and spacetime casual and geodesic
structure in modified gravity.

The energy conditions are the necessary conditions for a good understanding of
the singularity theorem, such as black-hole thermodynamics. The well known
Raychaudhuri equations (Raychaudhuri, 1955) play a role in describing the
attractive nature of gravity and positive energy density. The four basic
conditions are the null energy condition (NEC), weak energy condition (WEC),
dominant energy condition (DEC), and strong energy condition (SEC). All the
conditions are derived from the Raychaudhuri equation, which helps us analyze
the entire spacetime structure without precise solutions to Einstein’s
equations playing a vital role in understanding cosmological gravitational
interactions. The null energy condition (NEC) discusses the second law of
black-hole thermodynamics, although its violation corresponds to the
universe’s big rip singularity (Carroll, 2019). On the other strong energy
condition (SEC) is useful for studying the Hawking-Penrose theorem of
singularity (Hawking Ellis, 1973). SEC is also good at describing the
repulsive/attractive nature of gravity under modified theories of gravity. The
violation of SEC implies the observed accelerated expansion of the universe.
Many works on energy conditions under modified theories are presented in the
literature. K. Bamba et al. (Bamba et al, 2017) studied the energy conditions
in $f(G)$ gravity. They considered some realistic $f(G)$ forms that could be
used in the late-time cosmic acceleration epochs to cure finite-time future
singularities. S. Capozziello et al. present the role of energy conditions in
$f(R)$ cosmology (Capozziello ., 2018). K. Atazadeh et al. (Atazadeh Darabi,
2014) also considered two forms of $f(R,G)$ gravity accounting the stability
of the cosmological solutions using various energy conditions. M. Zubair et
al. (Zubair Waheed, 2015) worked on energy conditions in $f(T)$ gravity with
non-minimal torsion-matter coupling etc.

As it is known, modified gravity theories can be studied in two separate
groups based on a coupling between matter and geometry, i.e., minimal coupling
and non-minimal coupling. In general, the minimal coupling is weak, whereas
non-minimal coupling is a strong coupling. Several works investigate the
cosmological implications of a coupling of non-minimal matter-geometry in
$f(R,T)$. P.H.R.S. Moraes and P.K. Sahoo (Moraes Sahoo, 20172) studied the
most straightforward non-minimal matter-geometry coupling in the $f(R,T)$
cosmology. L.K. Sharma et al. (Sharma ., 2018) also worked on the existence of
non-minimal matter-geometry coupling govern by power-law in Bianchi I
spacetime. Recently, P.K. Sahoo and S. Bhattacharjee (Sahoo Bhattacharjee,
2020) investigated gravitational Baryogenesis in non-minimal coupled $f(R,T)$
gravity. Furthermore, most of the above studies are done by the prior
assumptions of scale factor or followed by the parametrization techniques.
However, working on these types of methods induce some unclear doubts on
constraining the model parameters for the current value of the universe’s age
and cosmological quantities. To overcome this issue, we have been adopted the
cosmography to study the present accelerated scenario of the universe in our
work.

In this work, we considered non-minimally coupled $f(R,T)$ gravity i.e.
$f(R,T)=R+\alpha RT$ where $\alpha$ is the model parameter, to study various
energy conditions. Also, to check the viability of $f(R,T)$ gravity theory, we
used the present values of the cosmological parameters $q_{0}$ and $H_{0}$.
The equation of state parameter $\omega$ provides an acceptable candidate for
comparing our models with $\Lambda$CDM. The recent results from Planck
Collaboration (Aghanim et al., 2020) and the $\Lambda$CDM model indicate that
the equation of state parameter $\omega\simeq-1$. This behavior corresponds to
the Universe’s negative pressure framework, which specifies the current
accelerated phase.

The literature includes several works with certain linear forms of $f(R,T)$
gravity (Sahoo ., 2020, 2016; Harko et al, 2011). We have considered the
linear form of $f(R,T)$ gravity in the last second section, which demonstrates
some consistency with $\Lambda$CDM and non-linear form. Using restricted
values of cosmological parameters, energy conditions are examined, which is
different from work done earlier.

The paper is presented and organized as follows: In section 2, we briefly
describe the formulation of field equations in $f(R,T)$ gravity. In the next
section, 3, we discuss the solutions of the field equations. The various
energy conditions are studied in section 4. Section 4 presents the comparison
with the $\Lambda$CDM model along with the equation of state parameter. Also,
the compatibility of model in case of linear form is presented in section 5.
Further discussions and conclusion are presented in section 6.

##  2 Basic equations of $f(R,T)$ gravity

The total action in the $f(R,T)$ theory of gravity reads (Harko et al, 2011)

$S=\frac{1}{16\pi}\int d^{4}x\sqrt{-g}f(R,T)+\int d^{4}x\sqrt{-g}L_{m},$ (1)

with $g$ being the metric determinant and $L_{m}$ be the matter Lagrangian
density, and the stress-energy tensor of matter is defined as

$T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}L_{m}\right)}{\delta
g^{\mu\nu}}.$ (2)

By varying this action (1) with respect to the metric yields

$\displaystyle\Big{[}f_{1}^{\prime}(R)$ $\displaystyle+$ $\displaystyle
f_{2}^{\prime}(R)f_{3}(T)\Bigr{]}R_{\mu\nu}-\frac{1}{2}f_{1}(R)g_{\mu\nu}$
$\displaystyle+$
$\displaystyle\Bigl{(}g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu}\Bigr{)}\Bigl{[}f_{1}^{\prime}(R)+f_{2}^{\prime}(R)f_{3}(T)\Bigr{]}$
$\displaystyle=\Bigl{[}8\pi$ $\displaystyle+$ $\displaystyle
f_{2}(R)f_{3}^{\prime}(T)\Bigr{]}T_{\mu\nu}+f_{2}(R)\left[f_{3}^{\prime}(T)p+\frac{1}{2}f_{3}(T)\right]g_{\mu\nu},$

for which it was assumed $f(R,T)=f_{1}(R)+f_{2}(R)f_{3}(T)$ and primes denote
derivatives with respect to the argument.

Now, we will take $f_{1}(R)=f_{2}(R)=R$ and $f_{3}(T)=\alpha T$, with $\alpha$
a constant. This is the simplest non-trivial functional form of the function
$f(R,T)$ which involves matter-geometry coupling within the $f(R,T)$
formalism. Moreover, it benefits from the fact that General Relativity is
retrieved when $\alpha=0$. Eq.(2) now reads

$G_{\mu\nu}=8\pi T_{\mu\nu},$ (4)

with $T_{\mu\nu}$ the usual matter energy-momentum tensor.

##  3 The $f(R,T)=R+\alpha RT$ cosmology

For a flat Friedmann-Robertson-Walker universe with scale factor $a(t)$ and
Hubble parameter $H=\dot{a}/a$, the expressions for energy density and
pressure from Eq. (2) as obtained in (Moraes Sahoo, 20172) reads as

$\rho\\!=\\!\frac{H^{2}\\!\left[8\pi-27\alpha\left(\dot{H}+2H^{2}\right)\right]\\!+\\!7\alpha(2\dot{H}+3H^{2})\left(\dot{H}+2H^{2}\right)}{\frac{64\pi^{2}}{3}-96\pi\alpha\left(\dot{H}+2H^{2}\right)+18\alpha^{2}\left(\dot{H}+2H^{2}\right)^{2}},$
(5)

$p\\!=\\!-\frac{9\alpha
H^{2}\\!\left(\dot{H}+2H^{2}\right)\\!+\\!\left(2\dot{H}\\!+\\!3H^{2}\right)\\!\left[\frac{8\pi}{3}-3\alpha\left(\dot{H}+2H^{2}\right)\right]}{\frac{64\pi^{2}}{3}-96\pi\alpha\left(\dot{H}+2H^{2}\right)+18\alpha^{2}\left(\dot{H}+2H^{2}\right)^{2}}.$
(6)

The study of cosmological scenarios is done by researchers using various
methods. Like the parametrization technique, this method is a good method to
display the universe’s profiles through mathematical modeling. Nevertheless,
some doubts about constraining the model parameters against the observational
data for present-day values of cosmological quantities are unclear. Keeping
this in mind, one can adopt the cosmography, which can overcome the previous
issues. The cosmographic parameters can be derived from the expansion of scale
factor concerning cosmic time. The Taylor series expansion of the scale factor
up to its fifth-order is given by

$\displaystyle a(t)$ $\displaystyle=$ $\displaystyle
a(t_{0})\Bigl{\\{}H_{0}(t-t_{0})-\frac{q_{0}}{2}H_{0}^{2}(t-t_{0})^{2}+\frac{j_{0}}{3!}H_{0}^{3}(t-t_{0})^{3}$
(7) $\displaystyle+$
$\displaystyle\frac{s_{0}}{4!}H_{0}^{4}(t-t_{0})^{4}+\frac{l_{0}}{5!}H_{0}^{5}(t-t_{0})^{5}+\mathcal{O}[(t-t_{0})^{6}]\Bigr{\\}},$

where the cosmological parameters such as Hubble parameter, deceleration
parameter, jerk, snap, and lerk can be written as

$\displaystyle
H=\frac{1}{a}\frac{da}{dt};\,\,\,\,\,\,q=-\frac{1}{aH^{2}}\frac{d^{2}a}{dt^{2}};$
$\displaystyle
j=\frac{1}{aH^{3}}\frac{d^{3}a}{dt^{3}};\,\,\,\,\,\,s=\frac{1}{aH^{4}}\frac{d^{4}a}{dt^{4}};$
$\displaystyle l=\frac{1}{aH^{5}}\frac{d^{5}a}{dt^{5}}.$ (8)

Also, we can customize the above terms as follows

$\displaystyle\dot{H}$ $\displaystyle=$ $\displaystyle-H^{2}(1+q)\,;$ (9)
$\displaystyle\ddot{H}$ $\displaystyle=$ $\displaystyle H^{3}(j+3q+2)\,;$ (10)
$\displaystyle\dddot{H}$ $\displaystyle=$ $\displaystyle
H^{4}[s-4j-3q(q+4)-6]\,.$ (11)

Now, using Equations (9), (10) in (5) and (6), we got the following
expressions

$\rho=\frac{3\alpha H^{4}(q-1)(7q+10)+12\pi
H^{2}}{27\alpha^{2}H^{4}(q-1)^{2}+144\pi\alpha H^{2}(q-1)+32\pi^{2}}\,;$ (12)
$p=\frac{H^{2}\left(9\alpha
H^{2}\left(q^{2}-1\right)+\pi(8q-4)\right)}{27\alpha^{2}H^{4}(q-1)^{2}+144\pi\alpha
H^{2}(q-1)+32\pi^{2}}\,.$ (13)

##  4 Energy Conditions

The energy conditions play a vital role in understanding the geodesics of the
Universe. These conditions are derived from the well-known Raychaudhuri
equations. For a congruence of time-like and null-like geodesics, the
Raychaudhuri equations are given in the following forms (Ehlers, 2006; Nojiri
Odinstov, 2007)

$\frac{d\theta}{d\tau}=-\frac{1}{3}\theta^{2}-\sigma_{\mu\nu}\sigma^{\mu\nu}+\omega_{\mu\nu}\omega^{\mu\nu}-R_{\mu\nu}u^{\mu}u^{\nu}\,,$
(14)

and

$\frac{d\theta}{d\tau}=-\frac{1}{2}\theta^{2}-\sigma_{\mu\nu}\sigma^{\mu\nu}+\omega_{\mu\nu}\omega^{\mu\nu}-R_{\mu\nu}n^{\mu}n^{\nu}\,,$
(15)

where $\theta$ is the expansion factor, $n^{\mu}$ is the null vector, and
$\sigma^{\mu\nu}$ and $\omega_{\mu\nu}$ are, respectively, the shear and the
rotation associated with the vector field $u^{\mu}$. For attractive gravity,
equations (18), and (19) satisfy the following conditions

$R_{\mu\nu}u^{\mu}u^{\nu}\geq 0\,;\,\,\,\,\,R_{\mu\nu}n^{\mu}n^{\nu}\geq 0\,.$
(16)

For perfect fluid matter distribution, the above conditions read the following

  * •

Strong energy conditions (SEC): $\rho+3p\geq 0$;

  * •

Weak energy conditions (WEC): $\rho\geq 0,\rho+p\geq 0$;

  * •

Null energy condition (NEC): $\rho+p\geq 0$;

  * •

Dominant energy conditions (DEC): $\rho\geq 0,|p|\leq\rho$.

Planck collaborations have recently improved the calculation of cosmological
parameters relative to the findings of 2015, leading to substantial
improvements in the accuracy of correlated parameters. The work assumed the
cosmology of the base $\Lambda$CDM and observed the parameters of the late
universe. So, it is found that the value of $H_{0}$ is $(67.4\pm 0.5)$ km s-1
Mpc-1 (Aghanim et al., 2020). Capozziello et al. also observed the
cosmographic parameters $q_{0},s_{0},j_{0}$ values (Capozziello ., 2019). The
cosmographic techniques are implemented to derive cosmological consistent
models of modified gravity. Further the value of $q_{0}$ is observed to be
$-0.503$.

To check the viability of $f(R,T)$ gravity theory we use the present values of
the cosmological parameters as $H_{0}=67.9$ and $q_{0}=-0.503$. The above
conditions reads as follows

$\rho+3p\\!=\\!\frac{3\alpha H_{0}^{4}(q_{0}-1)(16q_{0}+19)+24\pi
H_{0}^{2}q}{27\alpha^{2}H_{0}^{4}(q_{0}-1)^{2}+144\pi\alpha
H_{0}^{2}(q_{0}-1)+32\pi^{2}}\geq 0\,;$ (17) $\rho+p\\!=\\!\frac{3\alpha
H_{0}^{4}(q_{0}-1)(10q_{0}+13)+8\pi
H_{0}^{2}(q_{0}+1)}{27\alpha^{2}H_{0}^{4}(q_{0}-1)^{2}+144\pi\alpha
H_{0}^{2}(q_{0}-1)+32\pi^{2}}\geq 0\,;$ (18) $\rho\\!=\\!\frac{3\alpha
H_{0}^{4}(q_{0}-1)(7q_{0}+10)+12\pi
H_{0}^{2}}{27\alpha^{2}H_{0}^{4}(q_{0}-1)^{2}+144\pi\alpha
H_{0}^{2}(q_{0}-1)+32\pi^{2}}\geq 0\,;$ (19) $\rho-p\\!=\\!\frac{3\alpha
H_{0}^{4}(q_{0}-1)(4q_{0}+7)-8\pi
H_{0}^{2}(q_{0}-2)}{27\alpha^{2}H_{0}^{4}(q_{0}-1)^{2}+144\pi\alpha
H_{0}^{2}(q_{0}-1)+32\pi^{2}}\geq 0\,.$ (20)

Figure 1: Energy density $\rho$ for $f(R,T)=R+\alpha RT$ derived with the
present values of $H_{0}$ and $q_{0}$ parameters.

Figure 2: SEC $\rho+3p$ for $f(R,T)=R+\alpha RT$ derived with the present
values of $H_{0}$ and $q_{0}$ parameters.

Figure 3: NEC $\rho+p$ for $f(R,T)=R+\alpha RT$ derived with the present
values of $H_{0}$ and $q_{0}$ parameters.

Figure 4: DEC $\rho-p$ for $f(R,T)=R+\alpha RT$ derived with the present
values of $H_{0}$ and $q_{0}$ parameters.

From the above expressions (17)-(20), one can quickly observe that the ECs
depends on the model parameters $\alpha$. However, we can not choose the value
of $\alpha$ arbitrarily. If we do that may cause a violation of the current
accelerated scenario of the Universe. Keeping this thing in mind, we have
manipulated the profiles of the ECs in Figs. 1-4. As one can see, WEC, NEC,
DEC are satisfied from Figs. 1, 3 & 4 respectively while SEC violated with the
present value of the Hubble parameter ($H_{0}$) and deceleration parameter
($q_{0}$) in Fig 2. Also, this is an agreement with the current scenario of
the Universe.

The $\Lambda$CDM is a broadly accepted cosmological model that describes the
observations of the Universe. Several comments and surveys, such as Planck,
WAMP, and Dark Energy Survey (DES), tested it in the last two decades. A
specific mapping of $f(R,T)$ mimics the $\Lambda$CDM. So, for $\alpha=0$, our
model recover the $\Lambda$CDM, and the ECs for this are given by

  * •

SEC : $6H^{2}q\geq 0\,;$

  * •

NEC : $2H^{2}(1+q)\geq 0\,;$

  * •

WEC : $3H^{2}\geq 0$ and $2H^{2}(1+q)\geq 0\,;$

  * •

DEC : $3H^{2}\geq 0$ and $2H^{2}(1+q)\geq 0$ or, $-2H^{2}(-2+q)\geq 0\,.$

By applying the present value of $H_{0}$ and $q_{0}$, one can check that the
NEC, WEC, DEC satisfied whereas SEC violated. And these are the proper
behavior for a standard model, which shows the accelerated expansion of the
Universe. Moreover, our model is showing the same profiles for ECs, as
demonstrated by $\Lambda$CDM. Besides, the recent observation, such as the
Planck collaboration and the $\Lambda$CDM, confirms the equation of state
parameter $\omega$ takes its value as $\omega\simeq-1$ (Aghanim et al., 2020).
Also, this is an agreement for accelerated expansion. Therefore, $\omega$ is a
suitable candidate to compare our model with the $\Lambda$CDM model. The
expression of $\omega$ for our model is given by

$\omega=\frac{p}{\rho}=\frac{9\alpha
H^{2}\left(q^{2}-1\right)+\pi(8q-4)}{3\alpha H^{2}(q-1)(7q+10)+12\pi}\,.$ (21)

In Fig. 5, we have shown the profile of $\omega$ by taking the present value
of $H_{0}$ and $q_{0}$. One can observe that $\omega$ depends entirely on the
model parameter, and its value is very close to $-1$, which is in agreement
with the recent observations. Moreover, our model shows the accelerated
expansion as good as $\Lambda$CDM.

Figure 5: Equation of state parameter $\omega$ for $f(R,T)=R+\alpha RT$
derived with the present values of $H_{0}$ and $q_{0}$ parameters.

##  5 Compatibility with Linear $f(R,T)$ Case

In the present section, we study the linear case of $f(R,T)$ model i.e.,
$f(R,T)=R+2\gamma T$ (Harko et al, 2011), where $\gamma$ is a constant. Many
individuals operated in $f(R,T)$ gravity in the linear case, representing
energy conditions. But, by constrained values of cosmological parameters
$H_{0}$ and $q_{0}$, we have attempted to research the energy conditions. The
work includes the present values of $H_{0}=67.9$ and $q_{0}=-0.503$.

From (2), one can obtain the following equations

$3H^{2}=(1+3\gamma)\rho-\gamma p,$ (22) $2\dot{H}+3H^{2}=-(1+3\gamma)p+\gamma
p.$ (23)

Using (22) and (23), the energy density, pressure, and equation of state
parameter can be obtained as,

$\rho=\frac{-2\dot{H}\gamma+3(1+2\gamma)H^{2}}{(1+3\gamma)^{2}-\gamma^{2}},$
(24)
$p=-\frac{2(1+3\gamma)\dot{H}+3(1+2\gamma)H^{2}}{(1+3\gamma)^{2}-\gamma^{2}},$
(25)
$\omega=\frac{p}{\rho}=\frac{2H^{2}(3\gamma+1)(q+1)-3H^{2}(2\gamma+1)}{3H^{2}(2\gamma+1)+2H^{2}\gamma(q+1)}.$
(26)

The various energy conditions reads as,

$\rho+3p=\frac{2H_{0}^{2}[4\gamma+(10\gamma+3)q_{0}]}{8\gamma^{2}+6\gamma+1}\geq
0,$ (27) $\rho+p=\frac{2H_{0}^{2}(q_{0}+1)}{2\gamma+1}\geq 0,$ (28)
$\rho=\frac{2H_{0}^{2}(1+q_{0})\gamma+3(1+2\gamma)H_{0}^{2}}{(1+3\gamma)^{2}-\gamma^{2}}\geq
0,$ (29) $\rho-p=-\frac{2H_{0}^{2}(q_{0}-2)}{4\gamma+1}\geq 0.$ (30)

Figure 6: Energy conditions for $f(R,T)=R+2\gamma T$ derived with the present
values of $H_{0}$ and $q_{0}$ parameters. Figure 7: Equation of state
parameter $\omega$ for $f(R,T)=R+2\gamma T$ derived with the present values of
$H_{0}$ and $q_{0}$ parameters.

From the Fig. 6, we can observe that NEC, WEC, and DEC’s energy conditions
fulfill their conditions for $-0.25\leq\gamma\leq-0.15$ while SEC violates its
requirement. The SEC violation represents the accelerated expansion of the
universe. As we know, the EOS parameter is an acceptable candidate to compare
the model with the $\Lambda$CDM, so the $\gamma$ range is selected with the
EOS behavior observed.

The EOS parameter $\omega$ is approximately equal to -1 according to the
Planck observations and $\Lambda$CDM. The Fig. 7 also depicts the behavior of
$\omega$, taking the current values of $H_{0}$ and $q_{0}$. The values are
very similar to -1, which is in line with the $\Lambda$CDM. So, $f(R,T)$’s
linear and non-linear form is in good agreement with $\Lambda$CDM, resulting
in an accelerated universe expansion.

##  6 Conclusions

This paper built a cosmological paradigm from the simplest non-minimal matter-
geometry coupling in the gravitational theory of $f(R,T)$. We have considered
a well-motivated $f(R,T)$ gravity model as $f(R,T)=R+\alpha RT$, where
$\alpha$ is the model parameter.

The primary motivation for such a theory of gravity is linked to its
consistency with spacetime casual and geodesic structure discussed by various
energy conditions. This analysis derived the null, the weak, the dominant, and
strong energy conditions. In particular, the SEC plays an essential role in
describing gravity’s attractive or repulsive nature under the modified theory
of gravity. Also, the present values of cosmological parameters $H_{0}$ and
$q_{0}$ are used to check the viability of $f(R,T)$ gravity theory.

According to present values, NEC, WEC, and DEC are observed to satisfy the
conditions derived from the Raychaudhuri equations, whereas SEC is violated.
The energy conditions established $0.000045\leq\alpha\leq 0.000075$
constraints to describe an accelerated expansion of the Universe.We also study
the linear case of $f(R,T)$ model i.e., $f(R,T)=R+2\gamma T$, where $\gamma$
is a constant and observe that NEC, WEC, and DEC’s energy conditions fulfil
their conditions for $-0.25\leq\gamma\leq-0.15$ while SEC violates its
requirement. The SEC violation represents the accelerated expansion of the
universe. We also compared our energy constraints with those from the
$\Lambda$CDM model. If we talk about the $\Lambda$CDM gravity, all energy
conditions are satisfied except SEC. Furthermore, the equation of state
parameter, derived from our model, is consistent with a current negative
pressure period, showing values close to -1. This activity also confirms the
explanation of $\Lambda$CDM for dark energy and recent experimental
observations. So, $f(R,T)$’s linear and non-linear form is in good agreement
with $\Lambda$CDM, resulting in an accelerated universe expansion.

## Acknowledgements

PKS and SA acknowledges Council of Scientific & Industrial Research (CSIR),
New Delhi, India, for financial support to carry out the Research project No.
03(1454)/19/EMR-II Dt.02/08/2019. SM acknowledges Department of Science &
Technology (DST), Govt. of India, New Delhi, for awarding Junior Research
Fellowship file No. DST/INSPIRE Fellowship/2018/IF180676.

## References

  * Aghanim et al. (2020) PlanckAghanim, A. et al.  2020,  A & A641A6.
  * Alvarenga et al. (2013) Alvarenga2013Alvarenga, F. et al.  2013,  Phys. Rev. D87103526\.
  * Atazadeh Darabi (2014) Atazadeh2014Atazadeh, K. Darabi, F.  2014,  Gen. Relativ. Gravit.461664\.
  * Bamba et al (2017) Bamba2017Bamba, K. et al.  2017,  Gen. Relativ. Gravit.49112\.
  * Capozziello . (2019) Capozziello2019Capozziello, S., Agostino, R. Luongo, O.  2019,  Int. J. Mod. Phys. D281930016\.
  * Capozziello . (2018) Capo2018Capozziello, S., Nojiri, S. Odinstov, S.  2018,  Phys. Lett. B78199\.
  * Carroll (2019) Carroll2004Carroll, S.  2019,  Spacetime and Geometry: An Introduction to General Relativity Spacetime and Geometry: An Introduction to General Relativity.  CambridgeCambridge University Press.
  * Ehlers (2006) Ehlers2006Ehlers, J.  2006,  Int. J. Mod. Phys. D151573\.
  * Femaro Fiorini (2007) Femaro2007Femaro, R. Fiorini, F.  2007,  Phys. Rev. D75084031\.
  * Harko et al (2011) Harko2011Harko, T. et al.  2011,  Phys. Rev. D84024020\.
  * Hawking Ellis (1973) Hawking1973Hawking, S. Ellis, G.  1973,  The Large Scale Structure of Space-Time The Large Scale Structure of Space-Time.  CambridgeCambridge University Press.
  * Moraes Sahoo (20171) Moreas2017Moraes, P. Sahoo, P.  20171,  Phys. Rev. D96044038\.
  * Moraes Sahoo (20172) Pedro2017Moraes, P. Sahoo, P.  20172,  Eur. Phys. J. C77480\.
  * Nojiri Odinstov (2003) Nojiri2003Nojiri, S. Odinstov, S.  2003,  Phys. Rev. D68123512\.
  * Nojiri Odinstov (2007) Nojiri2007Nojiri, S. Odinstov, S.  2007,  Int. J. Geom. Methods Mod. Phys.4115\.
  * Nojiri . (2008) Nojiri2008Nojiri, S., Odinstov, S. Tretyakov, P.  2008,  Prog. Theor. Phys. Suppl.17281\.
  * Perlmutter et al. (1999) Perlmutter1999Perlmutter, S. et al.  1999,  Astron. J.1161009\.
  * Raychaudhuri (1955) Ray1955Raychaudhuri, A.  1955,  Phys. Rev.981123\.
  * Riess et al. (1998) Riess1998Riess, A. et al.  1998,  Astron. J.1161009\.
  * Riess et al. (1999) Riess1999Riess, A. et al.  1999,  Astron. J.117707\.
  * Sahoo Bhattacharjee (2020) Sahoo2020Sahoo, P. Bhattacharjee, S.  2020,  Int. J. Theor. Phys.591451\.
  * Sahoo . (2020) Psahoo2020Sahoo, P., Bhattacharjee, S., Tripathy, SK. Sahoo, PK.  2020,  Mod. Phys. Lett. A352050095\.
  * Sahoo . (2016) Sahoo2016Sahoo, P., Sahu, SK. Nath, A.  2016,  Eur. Phys. J. Plus13118\.
  * Sharif Zubair (2019) Sharif2012Sharif, M. Zubair, M.  2019,  J. Cosmol. Astropart. Phys.3028.
  * Sharma . (2018) Lokesh2018Sharma, L., Yadav, A., Sahoo, P. Singh, B.  2018,  Results in Physics10738\.
  * Yixin et al (2019) Yixin2019Yixin, X. et al.  2019,  Eur. Phys. J. C79708\.
  * Zubair Waheed (2015) Zubair2015Zubair, M. Waheed, S.  2015,  Astrophys. Space Sci.355361\.


# A projection method for porous media flow

Nicholas J. Derr  Chris H. Rycroft  John A. Paulson School of Engineering and
Applied Sciences, Harvard University, 29 Oxford Street, Cambridge, MA 02138
Mathematics Group, Lawrence Berkeley National Laboratory, 1 Cyclotron Road,
Berkeley, CA 94720

###### Abstract

Flow through porous, elastically deforming media is present in a variety of
natural contexts ranging from large-scale geophysics to cellular biology. In
the case of incompressible constituents, the porefluid pressure acts as a
Lagrange multiplier to satisfy the resulting constraint on fluid divergence.
The resulting system of equations is a possibly non-linear saddle-point
problem and difficult to solve numerically, requiring nonlinear implicit
solvers or flux-splitting methods. Here, we present a method for the
simulation of flow through porous media and its coupled elastic deformation.
The pore pressure field is calculated at each time step by correcting trial
velocities in a manner similar to Chorin projection methods. We demonstrate
the method’s second-order convergence in space and time and show its
application to phase separating neo-Hookean gels.

##  1 Introduction

Deforming porous media appear across nature, in settings such as soil
consolidation (Biot, 1941), river bed and channel formation (Rodríguez-Iturbe,
1997; Abrams et al., 2009), the planet’s mantle (McKenzie, 1984; Spiegelman,
1993), glacial melt runoff (Hewitt, 2011), living tissue (Ehlers et al., 2010;
Ranft et al., 2012), cellular interiors (Dembo, 1989; Howard, 2001; Bendix et
al., 2008; Schaller et al., 2013), and simple organisms (Radszuweit et al.,
2014). The essential similarity between all of these cases are the governing
behaviors laid out in the poroelasticity theory introduced by Biot (1941),
wherein a porous solid skeleton is suffused with fluid. Contraction of the
skeleton at a point in space raises the pore pressure, which gradually
equilibriates via fluid exchange with neighboring pores, leading to diffusion
of the pressure perturbation throughout the medium. Rice and Cleary (1976)
categorized the limiting behaviors of this response as undrained (on time
scales faster than fluid exchange) and drained (on times long enough for fluid
exchange to have taken place.) In contexts such as polymer gels, the fluid and
solid constituents are both effectively incompressible (Doi, 2009), such that
fluid exchange can be taken as effectively instantaneous and the system
behavior is effectively drained at all times.

Even in the linear regime, solving the coupled set of equations for pressure
and displacements is non-trivial. One widely used method for doing so is the
fixed-stress splitting method (Kim et al., 2011; Dana et al., 2022), where the
flow and mechanical problems are solved in succession. This is convenient in
that it allows separate code bases to be applied to individual parts of the
problem. However, this arrangement typically corresponds to obtaining
solutions on separate grids, necessitating coupling procedures and
interpolation between grids.

Many of the biological examples of porous media—especially polymer gels—are
capable of undergoing large volume changes that correspond to large
deformation and require non-linear constitutive models (Bendix et al., 2008;
Doi, 2009). In this regime, the equations of motion can be cast in either of
two different perspectives: Lagrangian, where variables are tracked at a given
material point; or Eulerian, where they are tracked at a given point in space.
Large deformation theories typically use a Lagrangian formulation for the
porestructure, tracking the displacement and fluid flux (relative to solid
motion) at each material point (De Boer, 2005). Alternatively, equivalent
Eulerian relations for the system can be derived from mixture theory as shown
by Coussy et al. (1998), tracking the density and momentum flux of each phase
at each point in space. These methods have previously been used to investigate
porous media flow as well (Spiegelman and McKenzie, 1987), relying on simple
geometries for analytical solutions.

In this paper, we propose a fully Eulerian method for porous media flow using
this latter formulation. While the method is general, we first intend to apply
it to the study of biological polymer gels. Individual monomers are
significantly less compressible than grains from rock or soil mechanics
(Howard, 2001), leading to the assumption above of constituent
incompressibility. This constraint is satisfied at all times by the pore
pressure field, acting as a Lagrange multiplier in much the same way as
pressure in the Navier–Stokes equations (De Boer and Didwania, 2001). This
corresponds to a non-linear multi-component saddle-point system which would
need to be solved several times per timestep, requiring specialized
preconditioners for the pressure block analogous to those required for
Brinkman flow (Borrvall and Petersson, 2003).

Instead, we satisfy the gel’s incompressibility constraint by adapting the
projection method for Navier–Stokes simulation, originally developed by Chorin
(1967, 1968). Here, we effectively apply a projection method (Bell et al.,
1989; Almgren et al., 1996; Puckett et al., 1997; Brown et al., 2001) to the
time-integration of the total mean velocity field obtained by volume averaging
those of each constituent. However, we do not use this mean velocity field as
a primary simulation variable as in the Lagrangian formulation above. At each
time step, an intermediate mean velocity field is calculated in the absence of
the constraint. This function is projected onto the space of incompressible
vector fields by solving a Poisson relation, and the result is used to correct
each constituent velocity field so that the incompressibility constraint is
satisfied.

In the next section, we review a derivation of the equations of porous media
flow, following the approach of Coussy et al. (1998). In section 3, we
describe the second-order spatial dicretizations used to approximate
derivatives and update rules for time-integration. In section 4, we apply the
method to two test cases: a manufactured solution for confirming the spatial
and temporal rates of convergence, and phase separation of viscoelastic
polymer gels.

##  2 Continuum theory

Consider a solid structure of characteristic pore size $\xi$ filled with
permeating fluid. We define continuum fields following the approach of Coussy
et al. (1998). On lengths comparable to the pore size—the microscale—each
point in space is occupied by polymer or solvent and has associated with it a
single value for fields such as density, momentum, or Cauchy stress state. On
larger lengths $L\gg\xi$—the macroscale—we can define corresponding fields at
a position $\mathbf{x}$ on the macroscale by averaging over a microscale
volume $V_{l}=\\{\mathbf{x}_{\mu}:|\mathbf{x}_{\mu}-\mathbf{x}|<l\\}$ centered
at $\mathbf{x}$ and of characteristic size $l$ with $L\gg l\gtrsim\xi$. In
particular, we can perform this averaging procedure for a single phase,
identified by an index $\alpha$, as

$\left<\psi\right>_{\alpha}\left(\mathbf{x},t\right)=\int_{V_{l}}\psi_{\mu}(\mathbf{x}_{\mu})K(\mathbf{x}_{\mu}-\mathbf{x})I_{\alpha}(\mathbf{x}_{\mu})\
dV.$ (1)

Here $\psi_{\mu}$ is a tensor or scalar field defined at the microscale, $K$
is a normalized kernel defined over balls of radius $l$, and $I_{\alpha}$ is
an indicator function equal to 1 if the point of space at $x_{i}^{\mu}$ is
occupied by constituent $\alpha$ and 0 otherwise. We consider for each phase
$\alpha$ a volume fraction $\phi_{\alpha}$, density $\rho_{\alpha}$, velocity
$\mathbf{v}_{\alpha}$, and Cauchy stress tensor $\bm{\uptau}_{\alpha}$ defined
according to

$\phi_{\alpha}=\left<1\right>_{\alpha},\qquad\rho_{\alpha}=\left<\rho_{\mu}\right>_{\alpha},\qquad\rho_{\alpha}\mathbf{v}_{\alpha}=\left<\rho_{\mu}\mathbf{v}_{\mu}\right>_{\alpha},\qquad\bm{\uptau}_{\alpha}=\left<\bm{\uptau}_{\mu}\right>_{\alpha}.$
(2)

Importantly, the densities defined here are apparent—corresponding to mass per
unit volume of total material rather than volume of constituent $\alpha$—and
the corresponding intrinsic densities can be written
$\hat{\rho}_{\alpha}=\rho_{\alpha}/\phi_{\alpha}$.

In the same way, we have defined apparent and intrinsic stress tensors. In
small-deformation regimes where $\phi$ is approximately constant, it is common
to close the problem by specifying the intrinsic stress
$\hat{\bm{\uptau}}_{\alpha}=\bm{\uptau}_{\alpha}/\phi_{\alpha}$ as a function
of the system variables (Radszuweit et al., 2014; Weber et al., 2018).
However, large deformation theories specify a form of the apparent stress
tensor based on the energetic cost of large-scale solid phase deformation (De
Boer, 2005). Derivatives and integrals can be defined in terms of similar
averaging procedures (Coussy et al., 1998). From this point forward, we will
concern ourselves only with relations at the macroscale.

In this work we will consider biphasic mixtures such that one phase is solid-
and the other fluid-like, in the sense of stresses depending on deformation
and deformation rate, respectively. Accordingly, we let
$\smash{\alpha=\text{f}}$ or s. Since $\phi_{\text{f}}+\phi_{\text{s}}=1$ by
definition, the volume fractions of both phases represent a single degree of
freedom; as such, we will frequently use $\phi=\phi_{\text{s}}$ to simply
denote the solid volume fraction.

###  2.1 Governing equations

At this point, we apply the principles of conservation of mass and
conservation of momentum to each phase in order to calculate the equations of
motion associated with gels (Doi, 2009). Conservation of mass within each
phase gives

$\frac{\partial\rho_{\alpha}}{\partial
t}+\bm{\nabla}\cdot\left(\rho_{\alpha}\mathbf{v}_{\alpha}\right)=s_{\alpha\leftarrow\beta}+s_{\alpha},$
(3)

where $s_{\alpha\leftarrow\beta}$ for $\alpha\neq\beta$ is the rate at which
mass of phase $\beta$ is transformed into mass of phase $\alpha$ and
$s_{\alpha}$ is the source rate of mass of phase $\alpha$. Conservation of
mass requires $s_{\alpha\leftarrow\beta}+s_{\beta\leftarrow\alpha}=0$. If the
constituents are incompressible such that the intrinsic density
$\hat{\rho}_{\alpha}$ is constant for each phase, this simplifies to a set of
volume fraction conservation equations,

$\frac{\partial\phi_{\alpha}}{\partial
t}+\bm{\nabla}\cdot\left(\phi_{\alpha}\mathbf{v}_{\alpha}\right)=\frac{s_{\alpha\leftarrow\beta}+s_{\alpha}}{\hat{\rho}_{\alpha}}.$
(4)

With substitution, (4) can be equivalently represented by a single volume
fraction conservation equation and an incompressibility condition on the mean
material velocity,

$\frac{\partial\phi}{\partial
t}+\bm{\nabla}\cdot\left(\phi\mathbf{v}_{\text{s}}\right)=\frac{s_{\text{s}\leftarrow\text{f}}}{\hat{\rho}^{\text{s}}},$
(5a)
$\bm{\nabla}\cdot\left[\phi\mathbf{v}_{\text{s}}+(1-\phi)\mathbf{v}_{\text{f}}\right]=\left(\frac{1}{\hat{\rho}^{\text{s}}}-\frac{1}{\hat{\rho}^{\text{f}}}\right)s_{\text{s}\leftarrow\text{f}}+\frac{s_{\text{s}}}{\hat{\rho}^{\text{s}}}+\frac{s_{\text{f}}}{\hat{\rho}^{\text{f}}}.$
(5b)

The first term on the right-hand side represents volume change due to phase
change (e.g. freezing or melting). The constraint (5b) necessitates the
introduction of a macroscopic variable $p$ to enforce it in the manner of a
Lagrange multiplier, exactly as in the case of the incompressible
Navier–Stokes equations. Going forward, we assume no material sources or phase
change, so that $s_{\text{s}\leftarrow\text{f}}=s_{\text{s}}=s_{\text{f}}=0$
and

$\bm{\nabla}\cdot\left[\phi\mathbf{v}_{\text{s}}+(1-\phi)\mathbf{v}_{\text{f}}\right]=0.$
(6)

Performing the same analysis using conservation of momentum rather than mass,
we obtain a momentum balance equation

$\rho_{\alpha}\left(\frac{\partial\mathbf{v}_{\alpha}}{\partial
t}+\mathbf{v}_{\alpha}\cdot\bm{\nabla}\mathbf{v}_{\alpha}\right)=\bm{\nabla}\cdot\bm{\uptau}_{\alpha}+\mathbf{f}_{\alpha\leftarrow\beta}+\mathbf{f}_{\alpha},$
(7)

where $\mathbf{f}_{\alpha}$ is a body force acting only on phase $\alpha$ and
$\mathbf{f}_{\alpha\leftarrow\beta}+\mathbf{f}_{\beta\leftarrow\alpha}=\mathbf{0}$
are a pair of equal-and-opposite forces exerted one each phase by the other.
These originate from integrals of surface tractions over phase boundaries on
the microscale. Assuming the existence of an isotropic stress $-p\mathbf{I}$
everywhere on the microscale, it can be shown (Coussy et al., 1998) that there
exists a pressure contribution to $\mathbf{f}_{\alpha\leftarrow\beta}$ equal
to $p\bm{\nabla}\phi_{\alpha}$. We assume the remainder of the interaction
force stems from the viscous stresses applied to solid surface by flow of the
penetrating fluid through the porestructure. Accordingly, we let

$\mathbf{f}_{\alpha\leftarrow\beta}=p\partial_{i}\phi_{\alpha}+\Gamma(\phi)\left(\mathbf{v}_{\beta}-\mathbf{v}_{\alpha}\right),$
(8)

where $\Gamma(\phi)\sim\mu/\xi^{2}$ is a prefactor scaling with pore-scale
viscous forces. We will let $\Gamma=\Gamma_{0}\phi^{2}/(1-\phi)$, $\Gamma_{0}$
a constant, for consistency with the Carman–Kozeny relation describing flow
through packed spheres (Carman, 1937; Kozeny, 1927).

To reflect the effect of letting $p$ exist in both phases, we decompose the
Cauchy stress $\bm{\uptau}_{\alpha}$. We define the effective stress
$\bm{\upsigma}_{\alpha}=\bm{\uptau}_{\alpha}+p\mathbf{I}$ as the portion of
macroscale stress not originating from pressure. We assume there are no other
sources of stress within the fluid phase besides the pressure $p$ and so set
$\bm{\upsigma}_{\text{f}}=\mathbf{0}$, i.e.
$\bm{\uptau}_{\text{f}}=-p\mathbf{I}$. This yields the identification of the
pressure $p$ as the porefluid pressure as in the theory of Biot (1941).
Finally, whatever the eventual definition of $\bm{\upsigma}_{\text{s}}$, we
assume the characteristic stress magnitudes are large relative to the inertial
forces on the left-hand side of (7), such that the left-hand side can be
omitted. In practice, we will preserve the partial time derivative and neglect
the advective term, reasoning that it is dominated by the right-hand side
terms.

With this notation and these assumptions, the set of model equations are

$\frac{\partial\phi}{\partial
t}+\bm{\nabla}\cdot(\phi\mathbf{v}_{\text{s}})=0,$ (9a)
$\rho_{\text{s}}\frac{\partial\mathbf{v}_{\text{s}}}{\partial
t}\approx\mathbf{0}=\bm{\nabla}\cdot\bm{\upsigma}_{\text{s}}-\phi\bm{\nabla}p+\Gamma(\phi)\left(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{s}}\right),$
(9b) $\rho_{\text{f}}\frac{\partial\mathbf{v}_{\text{f}}}{\partial
t}\approx\mathbf{0}=-(1-\phi)\bm{\nabla}p-\Gamma(\phi)\left(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{s}}\right),$
(9c)
$\bm{\nabla}\cdot\left[\phi\mathbf{v}_{\text{s}}+(1-\phi)\mathbf{v}_{\text{f}}\right]=0.$
(9d)

Note the two momentum balance equations can be summed to obtain a total
material mechanical equilibrium,

$\mathbf{0}=\bm{\nabla}\cdot\left(-p\mathbf{I}+\bm{\upsigma}_{\text{s}}\right)+\mathbf{f}_{\text{s}}+\mathbf{f}_{\text{f}},$
(10)

and in the absence of a body force $\mathbf{f}_{\text{f}}$, the fluid momentum
equation can be written in terms of the volumetric flux
$\mathbf{q}=(1-\phi)\left(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{s}}\right)$
as a Darcy-type relation

$\mathbf{q}=-\frac{(1-\phi)^{3}}{\Gamma_{0}\phi^{2}}\bm{\nabla}p,$ (11)

where we have substituted for $\Gamma(\phi)$, obtaining the usual
Carmen–Kozeny proportionality between flow rate and pressure gradients.

##  3 Numerical method

We now describe the numerical method by which we solve (9d). We begin by
discussing the discretization of our domain into a periodic, two-dimensional
grid. While we use the language of two dimensions to report the form of finite
difference stencils, we also consider their form in three dimensions to allow
for use of the method in either context.

###  3.1 Grid discretization

In two dimensions, we discretize a rectangular domain
$[a_{x},b_{x}]\times[a_{y},b_{y}]\in\mathbb{R}^{2}$ into an $m\times n$ grid
of cells $\Omega_{i,j}$ of width $\Delta x$ and height $\Delta y$. The left
bottom corner of cell $\Omega_{i,j}$ is

$\mathbf{c}_{i,j}=\left(a_{x}+i\Delta
x\right)\hat{\mathbf{x}}+\left(a_{y}+j\Delta y\right)\hat{\mathbf{y}},$ (12)

its center is

$\mathbf{x}_{i,j}=\mathbf{c}_{i,j}+\left(\frac{\Delta x\
\hat{\mathbf{x}}+\Delta y\ \hat{\mathbf{y}}}{2}\right),$ (13)

and the center of its left and bottom walls are given by

$\mathbf{w}_{i,j}^{l}=\mathbf{c}_{i,j}+\frac{\Delta y\
\hat{\mathbf{y}}}{2},\qquad\mathbf{w}_{i,j}^{d}=\mathbf{c}_{i,j}+\frac{\Delta
x\ \hat{\mathbf{x}}}{2}.$ (14)

We consider of series of time intervals $\Delta t$, such that we may write the
time of one step of the simulation as $t_{n}=n\Delta t$.

We discretize the simulation fields such that each cell is associated with a
single value of the solid fraction and phase velocities located at its center,

$\phi_{i,j}^{n}=\phi(\mathbf{x}_{i,j},t_{n}),\qquad(\mathbf{v}_{\text{s}})_{i,j}^{n}=\mathbf{v}_{\text{s}}(\mathbf{x}_{i,j},t_{n}),\qquad(\mathbf{v}_{\text{f}})_{i,j}^{n}=\mathbf{v}_{\text{f}}(\mathbf{x}_{i,j},t_{n}),$
(15)

Pressure and stress fields are located on cell corners and walls,
respectively, such that in each cell we have

$p_{i,j}^{n}=p(\mathbf{c}_{i,j},t_{n}),\qquad\bm{\upsigma}_{i,j}^{l,n}=\bm{\upsigma}(\mathbf{w}_{i,j}^{l},t_{n}),\qquad\bm{\upsigma}_{i,j}^{d,n}=\bm{\upsigma}(\mathbf{w}_{i,j}^{d},t_{n}).$
(16)

This allows for the simple calculation of finite difference approximations to
field derivatives. In particular, the pressure gradient at cell centers is
calculated

$\left(\bm{\nabla}p\right)^{n}_{i,j}=\left[\frac{p^{n}_{i+1,j+1}+p^{n}_{i+1,j}-p^{n}_{i,j+1}-p^{n}_{i,j}}{2\Delta
x}\right]\hat{\mathbf{x}}+\left[\frac{p^{n}_{i+1,j+1}+p^{n}_{i,j+1}-p^{n}_{i+1,j}-p^{n}_{i,j}}{2\Delta
y}\right]\hat{\mathbf{y}},$ (17)

and the stress divergence at cell centers is

$\left(\bm{\nabla}\cdot\bm{\upsigma}\right)^{n}_{i,j}=\frac{\hat{\mathbf{x}}\cdot\bm{\upsigma}^{l,n}_{i+1,j}-\hat{\mathbf{x}}\cdot\bm{\upsigma}^{l,n}_{i,j}}{\Delta
x}+\frac{\hat{\mathbf{y}}\cdot\bm{\upsigma}^{d,n}_{i,j+1}-\hat{\mathbf{y}}\cdot\bm{\upsigma}^{d,n}_{i,j}}{\Delta
y}.$ (18)

We also calculate the solid flux divergence at cell centers first by
introducing its value on the left and bottom walls,

$\left(\phi
v_{s}\right)^{l,n}_{i,j}=\frac{\left(\phi^{n}_{i-1,j}+\phi^{n}_{i,j}\right)\left(\hat{\mathbf{x}}\cdot(\mathbf{v}_{\text{s}})^{n}_{i-1,j}+\hat{\mathbf{x}}\cdot(\mathbf{v}_{\text{s}})^{n}_{i,j}\right)}{4},\\\
$ $\left(\phi
v_{s}\right)^{d,n}_{i,j}=\frac{\left(\phi^{n}_{i,j-1}+\phi^{n}_{i,j}\right)\left(\hat{\mathbf{y}}\cdot(\mathbf{v}_{\text{s}})^{n}_{i,j-1}+\hat{\mathbf{y}}\cdot(\mathbf{v}_{\text{s}})^{n}_{i,j}\right)}{4},$
(19a)

and then calculating the divergence analogously to the stress,

$\left(\bm{\nabla}\cdot(\phi\mathbf{v}_{\text{s}}\right))^{n}_{i,j}=\frac{\left(\phi
v_{s}\right)^{l,n}_{i+1,j}-\left(\phi v_{s}\right)^{l,n}_{i,j}}{\Delta
x}+\frac{\left(\phi v_{s}\right)^{d,n}_{i,j+1}-\left(\phi
v_{s}\right)^{d,n}_{i,j}}{\Delta y}.$ (20)

###  3.2 Arbitrary finite difference expressions

Now, we introduce second-order spatially converging finite difference
expressions for calculating discretized representations of various
differential operators. We will explicitly write out the finite difference
rules in two dimensions, letting $\psi_{i,j}=\psi(\mathbf{x}_{i,j})$ denote an
arbitrary cell-centered scalar field and
$\bm{\uppsi}_{i,j}=\bm{\uppsi}(\mathbf{x}_{i,j})$ an arbitrary vector field,
with $\bm{\uppsi}=\psi_{x}\hat{\mathbf{x}}+\psi_{y}\hat{\mathbf{y}}$.

First, we discuss the calculation of gradients of cell-centered fields on cell
walls. Writing the two-dimensional gradient out in terms of its components,
$\left(\bm{\nabla}\psi\right)^{w}=(\partial_{x}\psi)^{w}\hat{\mathbf{x}}+(\partial_{y}\psi)^{w}\hat{\mathbf{y}}$,
we will introduce separate finite difference rules for the normal derivatives
$(\partial_{x}\psi)^{l}$ and $(\partial_{y}\psi)^{d}$ and the tangential
derivatives $(\partial_{y}\psi)^{l}$ and $(\partial_{x}\psi)^{d}$. The normal
derivatives across the wall can be calculated by a simple difference as

$(\partial_{x}\psi)^{l}=\frac{\psi_{i,j}-\psi_{i-1,j}}{\Delta
x},\qquad(\partial_{y}\psi)^{d}=\frac{\psi_{i,j}-\psi_{i,j-1}}{\Delta y}.$
(21)

We compute tangential derivatives by averaging the differentiated quantity on
the cell walls above and below the gradient location in the derivative’s
direction, i.e.

$\left(\partial_{y}\psi\right)^{l}=\frac{\psi_{i,j+1}+\psi_{i-1,j+1}-\psi_{i,j-1}-\psi_{i-1,j-1}}{4\Delta
y},$ (22a)
$\left(\partial_{x}\psi\right)^{d}=\frac{\psi_{i+1,j}+\psi_{i+1,j-1}-\psi_{i-1,j}-\psi_{i-1,j-1}}{4\Delta
x}.$ (22b)

While we have written these expressions in terms of an arbitrary scalar field,
they also give the finite difference rules for the calculation of gradients of
vectors. In particular, for
$(\bm{\nabla}\bm{\uppsi})^{w}=\hat{\mathbf{x}}(\partial_{x}\bm{\uppsi})^{w}+\hat{\mathbf{y}}(\partial_{y}\bm{\uppsi})^{w}$,
the rules for the calculation of $(\partial_{x}\bm{\psi})^{w}$ and
$(\partial_{y}\bm{\psi})^{w}$ are exactly as above, with the vector field
$\bm{\uppsi}$ substituted for the scalar field $\psi$. This implies we
identify the row of a gradient tensor with the derivative index unless
otherwise specified.

The finite difference expressions can be generalized to three dimensions. In
two dimensions each wall features one normal and one tangential derivative,
but in three dimensions each wall features one normal and two tangential
derivatives.

Now we discuss the application of the advective derivative
$\mathbf{v}_{\alpha}\cdot\bm{\nabla}$ using an essentially non-oscillatory
(ENO) method which is stable in the face of shocks in the differentiated
fields. We use the second-order ENO2 method (Osher and Sethian, 1988; Osher
and Shu, 1991). In brief, the derivative in each direction is calculated from
a three-node stencil drawn from a set of five possible values. The three
selected nodes are chosen so that the approximated derivative is either
centered or upwinded and so that the nodes feature the smoothest variation of
the differentiated field in terms of the second derivative magnitude.

For simplicity, we consider a one-dimensional scalar field $\psi$ and scalar
velocity $v$ located at nodes spaced a distance $h$ apart in the
$x$-direction. The ENO2 approximation to the quantity $vd\psi/dx$ is
calculated as follows. At one node, we consider the set of five values
$S=\\{\psi_{-2},\psi_{-1},\psi_{0},\psi_{1},\psi_{2}\\}$ corresponding to the
field’s value at that cell and its first two neighbors in direction. With
these values, we can calculate the quantities

$s_{\text{b}}=|\psi_{-2}-2\psi_{-1}+\psi_{0}|,\qquad
s_{\text{c}}=|\psi_{-1}-2\psi_{0}+\psi_{1}|,\qquad
s_{\text{f}}=|\psi_{0}-2\psi_{1}+\psi_{2}|,$ (23)

which are proportional to the absolute magnitude of the second derivative at
the center cell as calculated using backward, centered, and forward second-
order stencils. We will use these as a proxy for which set of three nodes
feature sharper variation of the function $\psi$. We can also calculate three
approximations to the first derivative $f$ using each set of three nodes,

$f_{\text{b}}=\frac{\psi_{-2}-4\psi_{-1}+3\psi_{0}}{2h},\qquad
f_{\text{c}}=\frac{\psi_{1}-\psi_{-1}}{2h},\qquad
f_{\text{f}}=\frac{-3\psi_{0}+4\psi_{1}-\psi_{2}}{2h}.$ (24)

The approximation chosen for the first derivative is as follows: it uses
either the centered or upwinded stencil, choosing between the two based on
which has a lower sharpness as defined above. In other words, denoting the
ENO2 approximation to the first derivative as $E[S;v,h]$, we can write

$E\left[S;v,h\right]=\begin{cases}f_{\text{b}},&v>0\text{ and
}s_{\text{b}}<s_{\text{c}},\\\ f_{\text{f}},&v<0\text{ and
}s_{\text{f}}<s_{\text{c}},\\\ f_{\text{c}},&\text{otherwise.}\end{cases}$
(25)

Now, we explicitly relate this one-dimensional rule to the $D$-dimensional
vector operator applied to scalars and vectors in two dimensions. Consider a
scalar field discretized at cell centers such that
$\psi_{i,j}=\psi(\mathbf{x}_{i,j})$. If we denote
$v_{x}=\hat{\mathbf{x}}\cdot\mathbf{v}_{i,j}$,
$v_{y}=\hat{\mathbf{y}}\cdot\mathbf{v}_{i,j}$, then we can introduce the set
of five nodes in each dimension necessary for applying the ENO method,

$S^{(x)}_{i,j}=\\{\psi_{i-2,j},\psi_{i-1,j},\psi_{i,j},\psi_{i+1,j},\psi_{i+2,j}\\},\qquad
S^{(y)}_{i,j}=\\{\psi_{i,j-2},\psi_{i,j-1},\psi_{i,j},\psi_{i,j+1},\psi_{i,j+2}\\},$
(26)

and write the scalar advective derivative

$\left(\mathbf{v}\cdot\bm{\nabla}\psi\right)_{i,j}=v_{x}E\left[S^{(x)}_{i,j};v_{x},\Delta
x\right]+v_{y}E\left[S^{(y)}_{i,j};v_{y},\Delta y\right].$ (27)

The application to vector fields corresponds to independently applying the
operator to each scalar component field, i.e. for a vector field
$\bm{\uppsi}=\psi_{x}\hat{\mathbf{x}}+\psi_{y}\hat{\mathbf{y}}$, we write

$\left(\mathbf{v}\cdot\bm{\nabla}\bm{\uppsi}\right)_{i,j}=\left(\mathbf{v}\cdot\bm{\nabla}\psi_{x}\right)_{i,j}\hat{\mathbf{x}}+\left(\mathbf{v}\cdot\bm{\nabla}\psi_{y}\right)_{i,j}\hat{\mathbf{y}}.$
(28)

This differs from the gradient difference rule in that different stencils may
be chosen for the derivative of each component field, so it cannot be written
as a simple difference of vector values.

###  3.3 Update rules

Relations (9a)–(9c) each contain a partial time derivative, and as such are
amenable to explicit timestepping, but the incompressibility condition (9d)
contains no such update rule for the pressure. As previously mentioned, we
adapt Chorin-style projection methods to the system of equations presented
here (Chorin, 1967, 1968). We calculate intermediate velocity fields
$\mathbf{v}_{\text{s}}^{*}$ and $\mathbf{v}_{\text{f}}^{*}$ using the update
rules below that do not in general satisfy the incompressibility condition
(9d). Then, an elliptic PDE is solved to calculate the pressure required to
enforce incompressibility during the period $t\in[t_{n},t_{n+1}]$. As such,
the pressure obtained can be identified as that as the interval’s midpoint,
i.e.
$p^{\smash{n+\frac{1}{2}}}_{i,j}=p(\mathbf{c}_{ij},t_{n}+\frac{1}{2}\Delta
t)$, as detailed in e.g. by Brown et al. (2001). We include the previous time-
step’s pressure $p^{\smash{n-\frac{1}{2}}}_{i,j}$ in the calculation of the
intermediate velocities, so that the elliptic PDE correpsonds to the pressure
increment
$q^{n}_{ij}=p^{\smash{n+\frac{1}{2}}}_{ij}-p^{\smash{n-\frac{1}{2}}}_{ij}$.

We will decompose the stress tensor into elastic and viscous parts
$\bm{\upsigma}_{\text{s}}=\bm{\upsigma}_{e}+\bm{\upsigma}_{v}$. The portion of
the velocity update corresponding to the viscous stress will be completed with
the semi-implicit Crank–Nicolson method Crank and Nicolson (1947). Every other
aspect of the problem is explicitly integreated using the second-order Heun’s
method, a two-stage explicit Runge–Kutta method. We let $\overline{n+1}$
denote an approximate value at the next-time step computed with a forward
Euler update. An arbitrary field $\psi$ with $\partial\psi/\partial
t=g(\psi,t)$ would thus satisfy

$\psi^{\overline{n+1}}=\psi^{n}+\Delta tg(\psi,t_{n}),$ (29)

and a complete update with Heun’s method would give

$\psi^{n+1}=\psi^{n}+\frac{\Delta
t}{2}\left[g(\psi^{n},t_{n})+g(\psi^{\overline{n+1}},t_{n+1})\right].$ (30)

Thus, the solid fraction $\phi$ satisfies

$\frac{\phi^{n+1}_{i,j}-\phi^{n}_{i,j}}{\Delta
t}=-\frac{1}{2}\left[\left(\bm{\nabla}\cdot\phi\mathbf{v}_{\text{s}}\right)^{n}_{i,j}+\left(\bm{\nabla}\cdot\phi\mathbf{v}_{\text{s}}\right)^{\overline{n+1}}_{i,j}\right].$
(31)

Letting the drag force in cell $\Omega_{i,j}$ be denoted
$\bm{d}^{n}=\Gamma\left(\phi_{i,j}^{n})((\mathbf{v}_{\text{f}})_{i,j}^{n}-(\mathbf{v}_{\text{s}})_{i,j}^{n}\right)$,
we introduce solid and fluid velocity updates

$\frac{(\mathbf{v}_{\text{s}})^{*}_{i,j}-(\mathbf{v}_{\text{s}})^{n}_{i,j}}{\Delta
t}=\frac{1}{2}\left[\frac{\left(\bm{\nabla}\cdot\bm{\upsigma}_{e}\right)^{n}_{i,j}+\bm{d}^{n}}{\hat{\rho}_{\text{s}}\phi_{i,j}^{n}}+\frac{\left(\bm{\nabla}\cdot\bm{\upsigma}_{e}\right)^{\overline{n+1}}_{i,j}+\bm{d}^{\overline{n+1}}}{\hat{\rho}_{\text{s}}\phi_{i,j}^{\overline{n+1}}}\right]\\\
-\frac{\left(\bm{\nabla}p\right)^{n-\frac{1}{2}}_{i,j}}{\hat{\rho}_{s}}+\frac{1}{2}\left[\frac{\left(\bm{\nabla}\cdot\bm{\upsigma}_{v}\right)^{n}_{i,j}}{\hat{\rho}_{\text{s}}\phi^{n}_{i,j}}+\frac{\left(\bm{\nabla}\cdot\bm{\upsigma}_{v}\right)^{n+1}_{i,j}}{\hat{\rho}_{\text{s}}\phi^{n+1}_{i,j}}\right],$
(32)
$\frac{(\mathbf{v}_{\text{f}})^{*}_{i,j}-(\mathbf{v}_{\text{f}})^{n}_{i,j}}{\Delta
t}=\frac{1}{2}\left[\frac{-\bm{d}^{n}}{\hat{\rho}_{f}(1-\phi^{n}_{i,j})}+\frac{-\bm{d}^{\overline{n+1}}}{\hat{\rho}_{f}(1-\phi^{\overline{n+1}}_{i,j})}\right]-\frac{\left(\bm{\nabla}p\right)^{n-\frac{1}{2}}_{i,j}}{\hat{\rho}_{f}},$
(33)

where $\mathbf{v}_{\text{s}}^{*}$ and $\mathbf{v}_{\text{f}}^{*}$ are a set of
intermediate velocities that we do not expect to satisfy the condition of
incompressibility. From these we can compute an intermediate material velocity

$\bar{\mathbf{V}}^{*}_{i,j}=\phi^{*}_{i,j}(\mathbf{v}_{\text{s}})^{*}_{i,j}+(1-\phi^{*}_{i,j})(\mathbf{v}_{\text{f}})^{*}_{i,j}.$
(34)

At this point, we solve for the correction field $q^{n}_{i,j}$ such that
$p_{i,j}^{\smash{n+\frac{1}{2}}}=p^{\smash{n-\frac{1}{2}}}_{i,j}+q^{n}_{i,j}$.
For ease of notation, we introduce the specific volumes
$\nu_{\alpha}=1/\hat{\rho}_{\alpha}$. To indicate the following analysis is
without reference to any particular discretization, we omit the $(i,j)$
subscripts. We write the final velocity correction as

$\frac{\mathbf{v}_{\text{s}}^{n+1}-\mathbf{v}_{\text{s}}^{*}}{\Delta
t}=-\nu_{\text{s}}\bm{\nabla}q^{n},\qquad\frac{\mathbf{v}_{\text{f}}^{n+1}-\mathbf{v}_{\text{f}}^{*}}{\Delta
t}=-\nu_{\text{f}}\bm{\nabla}q^{n},$ (35)

The weighted average velocity $\bar{\mathbf{V}}=\phi\bm{u}+(1-\phi)\bm{v}$ at
timestep $n+1$ must be divergence-free. We construct the term
$\bm{\nabla}\cdot\bar{\mathbf{V}}^{n+1}$ by multiplying each of the above
relations by its corresponding phase fraction, applying the divergence
operator, and adding the two. Setting it to zero in accordance with the
incompressibility condition yields an elliptic PDE for the pressure increment
$q$,

$\bm{\nabla}\cdot\left[\Delta
t\bar{\nu}(\phi)\bm{\nabla}q^{n}\right]=\bm{\nabla}\cdot\bar{\mathbf{V}}^{*},$
(36)

where we have introduced
$\bar{\nu}(\phi)=\phi\nu_{\text{s}}+(1-\phi)\nu_{\text{f}}$, the material
specific volume.

Note the approximate $(n+1)$ velocity field in the Heun’s method time updates
correspond to the intermediate velocities, not the corrected divergence-free
velocities, introducing a local error $e^{n}\propto\Delta
t(\mathbf{v}^{n+1}_{\alpha}-\mathbf{v}^{*}_{\alpha})$ at each timestep.
However, this difference scales as
$(\mathbf{v}^{n+1}_{\alpha}-\mathbf{v}^{*}_{\alpha})\propto\Delta
t\bm{\nabla}q^{n}$, and the pressure increment $q^{n}$ scales with $\Delta t$
as well (Brown et al., 2001). Thus, $e^{n}\propto\Delta t^{3}$ and the second-
order convergence in time is preserved.

###  3.4 Numerical solution of the elliptic problem

Finally, we discuss the solution of (36) using the finite element method,
representing the pressure correction as a superposition of basis functions
$q(\mathbf{x},t_{n})=q_{i,j}^{n}\varphi_{i,j}(\mathbf{x})$. On each cell
$\Omega_{i,j}$, we introduce four basis functions. Letting
$c_{x}=\hat{\mathbf{x}}\cdot\mathbf{c}_{i,j}$ and
$c_{y}=\hat{\mathbf{y}}\cdot\mathbf{c}_{i,j}$, and letting
$\xi=(x-c_{x})/\Delta x$, $\eta=(y-c_{y})/\Delta y$,

$\varphi^{\text{dl}}_{i,j}=(1-\xi)(1-\eta),\
\varphi^{\text{dr}_{i,j}}=\xi(1-\eta),\
\varphi^{\text{ul}}_{i,j}=(1-\xi)\eta,\ \varphi^{\text{ur}}_{i,j}=\xi\eta,$
(37)

such that each basis function takes the value 1 at one corner of the cell and
0 at each of the others. We can then define the global basis function as

$\varphi_{i,j}(\mathbf{x})=\begin{cases}\varphi^{\text{dl}},&\mathbf{x}\in\Omega_{i,j},\\\
\varphi^{\text{dr}},&\mathbf{x}\in\Omega_{i-1,j},\\\
\varphi^{\text{ul}},&\mathbf{x}\in\Omega_{i,j-1},\\\
\varphi^{\text{ur}},&\mathbf{x}\in\Omega_{i-1,j-1},\\\
0,&\text{otherwise.}\end{cases}$ (38)

Multiplying by a test function $Q$ and integrating by parts, we obtain the
weak form of (36),

$\int_{V}\Delta t\bar{\nu}(\phi)\bm{\nabla}Q\cdot\bm{\nabla}q\
dV=\int_{V}\bar{\mathbf{V}}^{*}\cdot\bm{\nabla}Q\ dV-\int_{\partial
V}Q\hat{\mathbf{n}}\cdot\left(\bar{\mathbf{V}}^{*}-\Delta
t\bar{\nu}(\phi)\bm{\nabla}q\right)dS,$ (39)

but the surface integral on the right-hand side vanishes due to the
periodicity of our domain.

Letting a Greek index denote a single corner node of the grid, we write the
test function $Q=Q_{\nu}\varphi_{\nu}(\mathbf{x})$ and correction field
$q=q_{\nu}\varphi_{\nu}(\mathbf{x})$ in terms of our basis functions. The
above then corresponds to a linear system of equations determining the
$q$-coefficients,

$M_{\nu\gamma}q_{\gamma}=b_{\nu},$ (40)

with the operator $M$ and source vector $b$ defined as

$M_{\nu\gamma}=\int_{V}\Delta
t\bar{\nu}(\phi)\bm{\nabla}\varphi_{\alpha}\cdot\bm{\nabla}\varphi_{\beta}dV,\qquad
b_{\alpha}=\int_{V}\bar{\mathbf{V}}\cdot\bm{\nabla}\varphi_{\alpha}dV.$ (41)

At this point, we consider a model system in order to test the method. After
defining a stress tensor in the following section, we will be able to specify
a form for $\left(\bm{\nabla}\cdot\bm{\upsigma}\right)^{n}_{i,j}$ in (32),
which is the only undefined piece of information left to specify in the steps
above.

#### Stress specification

As mentioned, we decompose the stress tensor into elastic and viscous parts
$\bm{\upsigma}=\bm{\upsigma}_{e}+\bm{\upsigma}_{v}$. Elastic stresses are
produced by deformation of a solid body occupying the domain $\Omega$ relative
to an unstressed reference configuration $\Omega_{0}$. We will calculate
elastic stresses using the reference map technique (RMT) (Kamrin et al., 2012;
Rycroft et al., 2020; Lin et al., 2022). A position in the reference frame
$\bm{\xi}\in\Omega_{0}$ corresponds to a particular material point of the
deforming body, while a position in the deformed frame $\mathbf{x}\in\Omega$
corresponds to a particular location in physical space. The two are related by
a bijective mapping $\bm{\chi}[t]:\Omega_{0}\to\Omega$, such that the location
in physical space of material point $\bm{\xi}$ at time $t$ is
$\mathbf{x}(\bm{\xi},t)=\bm{\chi}[t](\bm{\xi})$. It is convenient in the case
of unstressed initial conditions to let the reference frame be the initial
body configuration, $\bm{\chi}[0](\bm{\xi})=\bm{\xi}$. We will use this
assumption going forward.

The deformation gradient tensor $\mathbf{F}$ is defined so that in indical
notation (and dropping for the moment the notion of time-dependence)
$F_{i\alpha}=\partial\chi_{i}/\partial\xi_{\alpha}$, where Latin and Greek
indices denote vector components in the deformed and reference frames,
respectively. Physically, $\mathbf{F}$ represents the mapping of material line
segments from the reference frame to the deformed frame, such that if
$\bm{d\chi}\in\Omega_{0}$ and $\bm{dx}\in\Omega$ refer to the same material
line segment in the two frames, they are related by
$dx_{i}=F_{i\alpha}d\chi_{\alpha}$. This quantity contains all of the
information relating to deformation of an elastic body. In particular, the
squared principal stretches $\lambda_{i}^{2}$ are given by the eigenvalues of
the left and right Cauchy-Green tensors
$\mathbf{B}=\mathbf{F}\cdot\mathbf{F}^{\mathsf{T}}$,
$\mathbf{C}=\mathbf{F}^{\mathsf{T}}\cdot\mathbf{F}$. In general we can write
$\bm{\upsigma}=\bm{\upsigma}(\mathbf{F})$ for a number of widely used finite-
strain models for elasticity.

To employ the RMT, we consider the reference map variable
$\mathbf{X}=\mathbf{X}(\mathbf{x},t)=\bm{\chi}[t]^{-1}(\mathbf{x})$ given by
the inverse of the bijective mapping $\bm{\chi}$. By tracking this field
relative to the constraint $D^{s}\mathbf{X}/Dt=0$, at each step and position
in space it is possible to calculate $\mathbf{F}^{-1}$ using our finite
difference stencil above, since $F^{-1}_{\alpha i}=\partial
X_{\alpha}/\partial x_{i}$. At the cost of an additional matrix inversion for
each calculation of the elastic stress tensor, we gain the ability to simulate
a wide variety of finite strain models for elastic stress on the same Eulerian
grid as the rest of our simulated fields. Here, we fix the simulated reference
map field $\mathbf{X}_{i,j}^{n}=\mathbf{X}(\mathbf{x}_{i,j},t_{n})$ at cell
centers and integrate the evolution equation $D^{s}\mathbf{X}/Dt=0$ using
Heun’s method,

$\frac{\mathbf{X}^{n+1}-\mathbf{X}^{n}}{\Delta
t}=-\frac{1}{2}\left[\left(\mathbf{v}_{\text{s}}\cdot\bm{\nabla}\mathbf{X}\right)^{n}_{i,j}+\left(\mathbf{v}_{\text{s}}\cdot\bm{\nabla}\mathbf{X}\right)^{\overline{n+1}}_{i,j}\right],$
(42)

where the advective derivatives are calculated using the ENO2 method. Given
the calculation of the deformation gradient as described, we let the solid
stress be given by a neo-Hookean model (Lin et al., 2022),

$\bm{\upsigma}_{e}=\kappa(J-1)\mathbf{I}+\frac{G}{J^{1+2/D}}\operatorname{dev}\left(\mathbf{F}\cdot\mathbf{F}^{\mathsf{T}}\right),\qquad
J=\operatorname{det}\mathbf{F},$ (43)

where $D$ is the dimension. Each of the quantities in the above can be
calculated from $\mathbf{F}$ as long as $\bm{\nabla}\mathbf{X}$ is evaluated
on cell walls using the finite difference rules established earlier. Note that
from $J=\operatorname{det}\mathbf{F}$ and (9a) we can derive

$J=\frac{\phi_{0}}{\phi},$ (44)

i.e. we can express the relative change of a differential volume element
either through the determinant of the deformation gradient $J$ or through the
inverse volume fraction $\phi$. Specifically, if the solid contracts ($J<0$)
the amount of solid at a position in space must go up ($\phi>\phi_{0}$) and
vice versa.

We assume there is also a dissipative stress related to the motion of the
polymer, as modeled by others including Radszuweit et al. (2014) and Weber et
al. (2018). Introducing the deformation rate tensor
$\mathbf{D}=\frac{1}{2}\left[\bm{\nabla}\mathbf{v}_{\text{s}}+(\bm{\nabla}\mathbf{v}_{\text{s}})^{\mathsf{T}}\right]$,
we write

$\bm{\upsigma}_{v}=\phi\eta_{b}\operatorname{tr}\left(\mathbf{D}\right)\mathbf{I}+2\phi\eta_{s}\operatorname{dev}\left(\mathbf{D}\right).$
(45)

In conjunction with the elastic stress $\bm{\upsigma}_{e}$, this yields a
Kelvin–Voigt viscoelastic material which behaves like a fluid on short times
and like a solid on long times. The Crank–Nicolson update in (32) is
constructed using the same finite difference stencils as in the explicit
updates.

It remains only to determine how we will solve the two linear systems arising
at each time step. For consistent implementation across two and three
dimensions, we solve each using PETSc, the Portable, Extensible Toolkit for
Scientific Computation (Balay et al., 1997, 2021a, 2021b). The rest of the
system is solved using custom-written C++ code with Open MPI (Gabriel et al.,
2004) for distributed-memory parallelization. The code base is templated on
dimension to minimize code propagation requirements and eliminate
implementation delay between the two- and three-dimensional cases.

##  4 Applications

In this section, we apply to the method above to two sample systems in order
to assess its performance, flexibility, and convergence.

###  4.1 Method of manufactured solutions

To verify the method’s rate convergence in space and time, we apply it to a
manufactured solution and report the error on several grid sizes. The
manufactured solution requires three parameters $A$, $P$ and $\bm{\uptau}$,
assumes $\hat{\rho}_{s}=\hat{\rho}_{f}=:\rho$, and is given by

$\phi(x,y,t)=\frac{1}{2}+A\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t},\qquad
p(x,y,t)=P\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t},$ (46)

where for ease of notation we have introduced the phase variables

$\varphi_{x}=\frac{2\pi x}{L_{x}},\ \varphi_{y}=\frac{2\pi y}{L_{y}},\
\varphi_{t}=\frac{2\pi t}{\bm{\uptau}}.$ (47)

This solid fraction is naturally produced by the velocity fields

$\mathbf{v}_{\text{s}}(x,y,t)=\left(\frac{-AL_{x}\sin\varphi_{x}\cos\varphi_{y}\cos\varphi_{t}}{\bm{\uptau}+2\bm{\uptau}A\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t}},\qquad\frac{-AL_{y}\cos\varphi_{x}\sin\varphi_{y}\cos\varphi_{t}}{\bm{\uptau}+2\bm{\uptau}A\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t}}\right)^{\mathsf{T}},$
(48)
$\mathbf{v}_{\text{f}}(x,y,t)=\left(\frac{AL_{x}\sin\varphi_{x}\cos\varphi_{y}\cos\varphi_{t}}{\bm{\uptau}-2\bm{\uptau}A\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t}},\qquad\frac{AL_{y}\cos\varphi_{x}\sin\varphi_{y}\cos\varphi_{t}}{\bm{\uptau}-2\bm{\uptau}A\cos\varphi_{x}\cos\varphi_{y}\sin\varphi_{t}}\right)^{\mathsf{T}}.$
(49)

Neglecting elasticity, the phase-specific body forces which will drive the
system towards the manufactures solution are

$\mathbf{f}_{s}(x,y,t)=\rho\phi\frac{\partial\mathbf{v}_{\text{s}}}{\partial
t}-\bm{\nabla}\cdot\bm{\upsigma}_{v}-\Gamma(\phi)(\mathbf{v}_{\text{f}}-\mathbf{v}_{\text{s}})+\phi\bm{\nabla}p,$
(50a)
$\mathbf{f}_{f}(x,y,t)=(1-\phi)\rho\frac{\partial\mathbf{v}_{\text{f}}}{\partial
t}-\Gamma(\phi)(\mathbf{v}_{\text{s}}-\mathbf{v}_{\text{f}})+(1-\phi)\bm{\nabla}p.$
(50b)

The resulting expressions are cumbersome, but a symbolic solver such as
Mathematica can be used to convert the expressions (50b) directly to C++ code.

(a) Manufactured solutionsolid fraction $\phi$ (b) Error scaling analysis
Figure 1: Scaling analysis using the method of manufactured solutions. (a) A
known solution, shown at $t=1.24$, manufactured using the body forces (50b)
with $A=0.1$, $P=0.1$, and $\bm{\uptau}=1$ over the domain
$\Omega=[0,10]^{2}$. The black lines indicate deformation of the solid phase
from $t=0$. (b) The system is simulated up to $t=10$, at which point the
numerical solution is compared with the actual values. The nodal $L^{2}$ error
is shown as a function of one-dimensional grid size, demonstrating second-
order convergence in space and time. The system parameters are
$\\{\rho,\eta_{b},\eta_{s},\Gamma_{0}\\}=\\{0.01,1,0.1,1\\}.$

A snapshot of the solution is shown on a $256^{2}$ grid in Figure 1(a). The
black lines in the figure corresponds to level sets of the reference map
components, and are plotted to illustrate deformation from the reference
configuration. To test the convergence, the system was simulated up to $t=10$
on square two-dimensional grids of varying sizes. The results are shown in
Figure 1(b), demonstrating second-order convergence in space and time.

###  4.2 Phase separating neo-Hookean gels

Polymer gels have been known for decades to phase separate under energetically
favorable conditions (Flory, 1953). Typically, this is described using a free
energy density of the form

$f_{\text{mix}}=\frac{k_{B}T}{\nu}\left[\phi\log\phi+(1-\phi)\log(1-\phi)+\chi\phi(1-\phi)\right],$
(51)

which represents the pointwise cost of maintaining a phase fraction $\phi$.
The logarithmic terms drive the system away from $\phi=0,1$ and are derived
from entropic considerations. The remaining term is scaled by the
Flory–Huggins interaction parameter, $\chi$, a scalar that represents the
polymer–solvent affinity or dis-affinity due to close-range energetic
interactions.

The total energy of the system is assumed to be $F=\int_{V}fdV$, with
$f=f_{\text{mix}}+\frac{k}{2}|\bm{\nabla}\phi|^{2}$, such that the presence of
boundaries between constant-$\phi$ phases is penalized. This is equivalent to
allowing for surface tension at the phase boundaries associated with an
interface width $\sqrt{k}$. Gradients in the chemical potential $\mu=\delta
f/\delta\phi$ drive flow of the solid phase relative to the fluid in order to
minimize $F$ over the system (Tanaka, 2000). This is represented as an
additional body force $\bm{f}_{t}$ acting on the solid phase given by

$\bm{f}_{t}=-\phi\bm{\nabla}\frac{\partial f}{\partial\phi}+\phi
k\bm{\nabla}\nabla^{2}\phi,$ (52)

which we incorporate into our system for the remainder of this section. We
simulate two 5122 systems with identical initial conditions
$\phi(x,y,0)=0.5+\delta\phi(x,y)$, where $\delta\phi(x,y)$ is a Gaussian
thermal noise field. The two systems differ only in the magnitude of the
elastic stress. Snapshots from each simulation are shown in Figure 2.

In Figure 2(a), the solid phase experiences no elastic stress, and we recover
the behavior of liquid–liquid phase separation, albeit with asymmetric
rheologies as in the work of Tanaka (1993). The initial instability proceeds
unimpeded, and the gel forms biconnected regions before settling into an array
of droplets which undergo slow ripening. In Figure 2(b), the elastic
resistance slows the rate at which phase separation occurs, and eventually the
dynamics are arrested as thermodynamic forces are exactly balanced by the
elastic reponse to the extreme deformation.

solid fraction $\phi$$t=180$$t=250$$t=500$$t=1000$$K=0.5$, $G=0.05$$K=G=0$
Figure 2: Comparison of classical phase separation with and without elastic
resistance. (a) In the absence of elasticity, a biphasic gel undergoes
spinodal decomposition, forming biconnected regions and eventually isolated
droplets which undergo slow ripening (Tanaka, 1993). (b) In the presence of
elasticity, phase separation proceeds on a slower time scale and is ultimately
arrested. As in Figure 1, the black lines in the second row are level sets of
the reference map field and show deformation of the solid phase. Note that
while in Figure 1 the lines are purely for visualization, in this case the
reference map field is used to calculate the elastic resistance opposing phase
separation. The other system parameters are
$\\{\rho,k_{b}T/\nu,\eta_{b},\eta_{s},\Gamma_{0},\chi,k\\}=\\{0.01,1,1,0.1,1,2.4,0.0625\\}.$

For this system, we also report how the method scales across differing numbers
of processors in two and three dimensions. Figure 3 shows the parallel
efficiency for systems of varying grid size across 1–16 CPU cores. For a given
system size, we record the serial wall clock time $T_{1}$ to simulate a
certain number of time steps $n$. We do not record time for initialization or
data recording. For a given number of cores $k$, we then compute the time
$T_{k}$ to simulate $n$ timesteps and calculate the efficiency as
$T_{1}/kT_{k}$. In two dimensions, the efficiency is not fully realized until
grid sizes around 2562 or above. In that regime, it reaches around 75%
efficiency on 2 cores and decreases to 60% on 16 cores. The three dimensional
case is generally more efficient, which is to be expected given the large
system sizes. At 2 cores, efficiencies are at or above 80%, and at 16 cores
they remain at 70–80%.

Figure 3:  Parallel efficiencies of the computational method on 1–16 CPU cores
for a variety of grid sizes in two (left) and three (right) dimensions. The
duration of simulation is constant for a given grid size. The efficiency
$e_{k}$ on $k$ cores is computed by comparing the wall time to simulate on a
single core $T_{1}$ and on $k$ cores $T_{k}$. Then $e_{k}=T_{1}/kT_{k}$. The
domain is $[0,10]$ in each dimension, and the system parameters are
$\\{\rho,\kappa,G,k_{b}T/\nu,\eta_{b},\eta_{s},\Gamma_{0},\chi,k\\}=\\{0.01,0,0,1,1,0.1,1,2.4,0.0625\\}$.

##  5 Conclusion

In this work, we have introduced a novel numerical method for the simulation
of incompressible gels which will be applied to the study of active
cytoskeletal gels. The method is completely Eulerian, leading to
straightforward coupling between flow and mechanics. Moreover, with this
approach it is easy to introduce additional model fields coupled to the solid
or fluid phase, allowing for the development of multi-physics methods. It is
general with respect to the solid stress and can thus be applied to linear,
non-linear, elastic, and visco-elastic porestructure rheologies.

We have shown the method converges at second-order in time and space, and by
using a semi-implicit update for the solid velocity phase we eliminate the
second-order CFL condition which would otherwise arise due to solid
viscoelastic behavior. By applying a projection step to enforce the material
incompressibility condition, we avoid difficult-to-solve saddle-point
structures. Similarly, we by using explicit updates and the reference map
method, we allow for simulation of non-linear elastic and viscoelastic
rheologies without requiring Newton solves at each time step.

Many porous media problems of interest feature boundaries, and the method is
amenable to specialization towards that end through the use of level sets or
the reference map itself as in (Lin et al., 2022). While the equations it
simulates are often derived in the context of soil and rock mechanics, the
method is well suited for describing porous systems across nature, including
biological settings such as developing tissue or cellular interiors.

## Acknowledgements

We thank Christoph Weber (Max Planck Institute for the Physics of Complex
Systems) and L. Mahadevan (Harvard University) for useful discussions about
this work. This work was supported by the National Science Foundation under
Grant No. DMS-1753203. C. H. R. was partially supported by the Applied
Mathematics Program of the U.S. DOE Office of Science Advanced Scientific
Computing Research under Contract No. DE-AC02-05CH11231. N. J. D. was
partially supported by the NSF–Simons Center for Mathematical and Statistical
Analysis of Biology at Harvard (Award No. 1764269) and the Harvard
Quantitative Biology Initiative.

## References

  * Biot (1941) M. A. Biot,  General theory of three‐dimensional consolidation,  Journal of Applied Physics 12 (1941) 155–164.
  * Rodríguez-Iturbe (1997) I. Rodríguez-Iturbe, Fractal River Basins: Change and Self-Organization, Cambridge University Press, 1997.
  * Abrams et al. (2009) D. M. Abrams, A. E. Lobkovsky, A. P. Petroff, K. M. Straub, B. McElroy, D. C. Mohrig, A. Kudrolli, D. H. Rothman,  Growth laws for channel networks incised by groundwater flow,  Nature Geoscience 2 (2009) 193–196.
  * McKenzie (1984) D. McKenzie,  The Generation and Compaction of Partially Molten Rock,  Journal of Petrology 25 (1984) 713–765.
  * Spiegelman (1993) M. Spiegelman,  Flow in deformable porous media. Part 1 Simple analysis,  Journal of Fluid Mechanics 247 (1993) 17–38.
  * Hewitt (2011) I. J. Hewitt,  Modelling distributed and channelized subglacial drainage: the spacing of channels,  Journal of Glaciology 57 (2011) 302–314.
  * Ehlers et al. (2010) W. Ehlers, A. Acartürk, N. Karajan,  Advances in modelling saturated soft biological tissues and chemically active gels,  Archive of Applied Mechanics 80 (2010) 467–478.
  * Ranft et al. (2012) J. Ranft, J. Prost, F. Jülicher, J. F. Joanny,  Tissue dynamics with permeation,  The European Physical Journal E 35 (2012) 46.
  * Dembo (1989) M. Dembo,  Mechanics and control of the cytoskeleton in Amoeba proteus,  Biophysical Journal 55 (1989) 1053–1080.
  * Howard (2001) J. Howard, Mechanics of Motor Proteins and the Cytoskeleton, Sinauer Associates, 2001\.
  * Bendix et al. (2008) P. M. Bendix, G. H. Koenderink, D. Cuvelier, Z. Dogic, B. N. Koeleman, W. M. Brieher, C. M. Field, L. Mahadevan, D. A. Weitz,  A Quantitative Analysis of Contractility in Active Cytoskeletal Protein Networks,  Biophysical Journal 94 (2008) 3126–3136.
  * Schaller et al. (2013) V. Schaller, K. M. Schmoller, E. Karaköse, B. Hammerich, M. Maier, A. R. Bausch,  Crosslinking proteins modulate the self-organization of driven systems,  Soft Matter 9 (2013) 7229.
  * Radszuweit et al. (2014) M. Radszuweit, H. Engel, M. Bär,  An Active Poroelastic Model for Mechanochemical Patterns in Protoplasmic Droplets of Physarum polycephalum,  PLoS ONE 9 (2014) e99220.
  * Rice and Cleary (1976) J. R. Rice, M. P. Cleary,  Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents,  Reviews of Geophysics 14 (1976) 227–241.
  * Doi (2009) M. Doi,  Gel Dynamics,  Journal of the Physical Society of Japan 78 (2009) 052001.
  * Kim et al. (2011) J. Kim, H. A. Tchelepi, R. Juanes,  Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits,  Computer Methods in Applied Mechanics and Engineering 200 (2011) 1591–1606.
  * Dana et al. (2022) S. Dana, M. Jammoul, M. F. Wheeler,  Performance studies of the fixed stress split algorithm for immiscible two-phase flow coupled with linear poromechanics,  Computational Geosciences 26 (2022) 13–27.
  * De Boer (2005) R. De Boer, Trends in Continuum Mechanics of Porous Media, volume 18, Springer, 2005\.
  * Coussy et al. (1998) O. Coussy, L. Dormieux, E. Detournay,  From mixture theory to Biot’s approach for porous media,  International Journal of Solids and Structures 35 (1998) 4619–4635.
  * Spiegelman and McKenzie (1987) M. Spiegelman, D. McKenzie,  Simple 2-D models for melt extraction at mid-ocean ridges and island arcs,  Earth and Planetary Science Letters 83 (1987) 137–152.
  * De Boer and Didwania (2001) R. De Boer, A. K. Didwania,  Saturated elastic porous solids: Incompressible, compressible and hybrid binary models,  Transport in Porous Media 45 (2001) 423–443.
  * Borrvall and Petersson (2003) T. Borrvall, J. Petersson,  Topology optimization of fluids in Stokes flow,  International Journal for Numerical Methods in Fluids 41 (2003) 77–107.
  * Chorin (1967) A. J. Chorin,  A numerical method for solving incompressible viscous flow problems,  Journal of Computational Physics 2 (1967) 12–26.
  * Chorin (1968) A. J. Chorin,  Numerical solution of the Navier–Stokes equations,  Mathematics of Computation 22 (1968) 745–762.
  * Bell et al. (1989) J. B. Bell, P. Colella, H. M. Glaz,  A second-order projection method for the incompressible Navier–Stokes equations,  Journal of Computational Physics 85 (1989) 257–283.
  * Almgren et al. (1996) A. S. Almgren, J. B. Bell, W. G. Szymczak,  A numerical method for the incompressible Navier–Stokes equations based on an approximate projection,  SIAM Journal on Scientific Computing 17 (1996) 358–369.
  * Puckett et al. (1997) E. G. Puckett, A. S. Almgren, J. B. Bell, D. L. Marcus, W. J. Rider,  A high-order projection method for tracking fluid interfaces in variable density incompressible flows,  Journal of Computational Physics 130 (1997) 269–282.
  * Brown et al. (2001) D. L. Brown, R. Cortez, M. L. Minion,  Accurate projection methods for the incompressible Navier–Stokes equations,  Journal of Computational Physics 168 (2001) 464–499.
  * Weber et al. (2018) C. A. Weber, C. H. Rycroft, L. Mahadevan,  Differential Activity-Driven Instabilities in Biphasic Active Matter,  Physical Review Letters 120 (2018) 248003.
  * Carman (1937) P. C. Carman,  Fluid flow through granular beds,  Transactions of the Institution of Chemical Engineers 15 (1937) 150–166.
  * Kozeny (1927) J. Kozeny,  Uber kapillare leitung der wasser in boden,  Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften. Mathematisch-Naturwissenschaftliche Classe 136 (1927) 271–306.
  * Osher and Sethian (1988) S. Osher, J. A. Sethian,  Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton–Jacobi formulations,  Journal of Computational Physics 79 (1988) 12–49.
  * Osher and Shu (1991) S. Osher, C.-W. Shu,  High-order essentially nonoscillatory schemes for Hamilton–Jacobi equations,  SIAM Journal on Numerical Analysis 28 (1991) 907–922.
  * Crank and Nicolson (1947) J. Crank, P. Nicolson,  A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,  Mathematical Proceedings of the Cambridge Philosophical Society 43 (1947) 50–67.
  * Kamrin et al. (2012) K. Kamrin, C. H. Rycroft, J.-C. Nave,  Reference map technique for finite-strain elasticity and fluid–solid interaction,  Journal of the Mechanics and Physics of Solids 60 (2012) 1952–1969.
  * Rycroft et al. (2020) C. H. Rycroft, C.-H. Wu, Y. Yu, K. Kamrin,  Reference map technique for incompressible fluid–structure interaction,  Journal of Fluid Mechanics 898 (2020) A9.
  * Lin et al. (2022) Y. L. Lin, N. J. Derr, C. H. Rycroft,  Eulerian simulation of complex suspensions and biolocomotion in three dimensions,  Proceedings of the National Academy of Sciences 119 (2022) e2105338118.
  * Balay et al. (1997) S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith,  Efficient management of parallelism in object oriented numerical software libraries,  in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
  * Balay et al. (2021a) S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, J. Zhang, PETSc/TAO Users Manual, Technical Report ANL-21/39 - Revision 3.16, Argonne National Laboratory, 2021a.
  * Balay et al. (2021b) S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, J. Zhang, PETSc Web page, https://petsc.org/, 2021b. URL: https://petsc.org/.
  * Gabriel et al. (2004) E. Gabriel, G. E. Fagg, G. Bosilca, T. Angskun, J. J. Dongarra, J. M. Squyres, V. Sahay, P. Kambadur, B. Barrett, A. Lumsdaine, et al.,  Open MPI: Goals, concept, and design of a next generation MPI implementation,  in: European Parallel Virtual Machine/Message Passing Interface Users’ Group Meeting, Springer, 2004, pp. 97–104.
  * Flory (1953) P. J. Flory, Principles of Polymer Chemistry, Cornell University Press, 1953\.
  * Tanaka (2000) H. Tanaka,  Viscoelastic phase separation,  Journal of Physics: Condensed Matter 12 (2000) R207–R264.
  * Tanaka (1993) H. Tanaka,  Unusual phase separation in a polymer solution caused by asymmetric molecular dynamics,  Physical Review Letters 71 (1993) 3158–3161.