$$\begin{equation} \begin{aligned} & D\delta W_{int}\left[\Delta\mathbf{u}\right]+D\delta W_{int}\left[\Delta\mathbf{w}\right]+D\delta W_{int}\left[\Delta J^{f}\right]\\ & -D\delta W_{ext}\left[\Delta\mathbf{u}\right]-D\delta W_{ext}\left[\Delta\mathbf{w}\right]-D\delta W_{ext}\left[\Delta J^{f}\right]\\ & \approx\delta W_{ext}-\delta W_{int}\,. \end{aligned} \label{eq:linearized-work-split-bfsi} \end{equation}$$

The linearizations of integrals are performed in the material frame of the solid domain of $\Omega^{b}$ , allowing us to linearize $\delta W_{int}$ along increments $\Delta\mathbf{u}$ , $\Delta\mathbf{w}$ , or $\Delta J^{f}$ by simply bringing the directional derivative operator inside the integrals of Eqs. (3.7.1-3)-(3.7.1-4). For notational convenience, we let $J\equiv J^{s}$ and $\mathbf{F}\equiv\mathbf{F}^{s}$ . Thus, the conversion of the internal virtual work to the material frame of the solid produces $$\begin{equation} \begin{aligned}\delta W_{int}= & \int_{\Omega^{b}}\left(-\varphi_{r}^{s}\mathbf{T}_{d}+\mathbf{F}\cdot\boldsymbol{\sigma}^{e}\cdot\mathbf{F}^{T}\right):\Grad\delta\mathbf{v}\cdot\mathbf{F}^{-1}dV\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\varphi_{r}^{s}\left(\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}\cdot\mathbf{F}^{-T}\cdot\Grad\left(\frac{\varphi^{f}}{\varphi^{s}}\right)-\rho_{T}^{f}\mathbf{a}^{f}+\rho_{T}^{s}\dot{\mathbf{v}}^{s}\right)dV\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-J^{2}\mathbf{F}^{-T}\cdot\mathbf{K}^{-1}\cdot\mathbf{F}^{-1}\cdot\mathbf{w}dV\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot J\left(\rho_{T}^{f}\mathbf{a}^{f}+\mathbf{F}^{-T}\cdot\Grad p\right)dV\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot J\left(-\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\mathbf{F}^{-T}\cdot\Grad\varphi^{f}+J\mathbf{F}^{-T}\cdot\mathbf{K}^{-1}\cdot\mathbf{F}^{-1}\cdot\mathbf{w}\right)dV\\ & +\int_{\Omega^{b}}J\boldsymbol{\tau}:\Grad\delta\mathbf{w}\cdot\mathbf{F}^{-1}dV\\ & +\int_{\Omega^{b}}J\Grad\delta J^{f}\cdot\mathbf{F}^{-1}\cdot\mathbf{w}+\delta J^{f}\left(J\frac{\varphi^{f}}{J^{f}}\frac{D^{f}J^{f}}{Dt}-\dot{J}^{s}\right)dV, \end{aligned} \label{eq:int-virtual-work-material-BFSI} \end{equation}$$ where $\mathbf{S}=J\,\mathbf{F}^{-1}\cdot\boldsymbol{\sigma}^{e}\cdot\mathbf{F}^{-T}$ is the second Piola-Kirchhoff stress for the solid constituent of the mixture, $\mathbf{W}=J\mathbf{F}^{-1}\cdot\mathbf{w}$ is the Piola transformation of $\mathbf{w}$ , and $dV=J^{-1}dv$ is an elemental volume of $\Omega^{f}$ in its material frame [23]. In addition, $\mathbf{K}^{-1}=J^{-1}\mathbf{F}^{T}\cdot\mathbf{k}^{-1}\cdot\mathbf{F}$ is the inverse of the permeability tensor in the material frame. Note that in the material frame, the fluid acceleration $\mathbf{a}^{f}$ is $$\begin{equation} \mathbf{a}^{f}=\frac{1}{\varphi^{f}}\left(\dot{\mathbf{w}}+\varphi^{f}\dot{\mathbf{v}}^{s}-\frac{\varphi^{s}}{\varphi^{f}}\frac{\dot{J}}{J}\mathbf{w}+\frac{1}{J}\left(\frac{1}{\varphi^{f}}\Grad\mathbf{w}+\Grad\mathbf{v}^{s}-\frac{1}{\varphi^{f^{2}}}\mathbf{w}\otimes\Grad\varphi^{f}\right)\cdot\mathbf{W}\right)\,.\label{eq:material-fluid-accel-BFSI} \end{equation}$$ The linearization of $\delta W_{int}$ is then performed along an increment $\Delta\mathbf{u}$ , and the integral is reverted back to the spatial frame, yielding for the displacement equations $$\begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{u}\right]= & \int_{\Omega^{b}}\alpha_{f}\left(\grad\Delta\mathbf{u}\cdot\boldsymbol{\sigma}^{e}:\grad\delta\mathbf{v}+\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}:\grad\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\varphi^{s}\boldsymbol{\tau}:\grad\delta\mathbf{v}\cdot\left(\grad\Delta\mathbf{u}+\frac{\varphi^{s}}{\varphi^{f}}\divg\Delta\mathbf{u}\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}^{\tau}:\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{D}^{w}+\mathbf{M}\cdot\grad\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\grad\delta\mathbf{v}:\frac{1}{\varphi^{f^{2}}}\boldsymbol{\mathcal{C}}^{\tau}:2\frac{\varphi^{s}}{\varphi^{f}}\left(\divg\Delta\mathbf{u}\right)\mathbf{w}\otimes\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\alpha_{f}\varphi^{s}\grad\delta\mathbf{v}:\frac{1}{\varphi^{f^{2}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\grad^{T}\Delta\mathbf{u}\cdot\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\varphi^{s}\grad\delta\mathbf{v}:\frac{1}{\varphi^{f^{2}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\left(\divg\Delta\mathbf{u}\right)\grad\varphi^{f}+\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\tau}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\divg\Delta\mathbf{u}+\grad^{T}\Delta\mathbf{u}\right)\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\tau}\cdot\frac{1}{\varphi^{s}}\grad\left(\divg\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\mathcal{C}}^{\tau}:\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{D}^{w}+\mathbf{M}\cdot\grad\Delta\mathbf{u}\right)\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:2\frac{\varphi^{s}}{\varphi^{f}}\left(\divg\Delta\mathbf{u}\right)\mathbf{w}\otimes\grad\varphi^{f}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\grad^{T}\Delta\mathbf{u}\cdot\grad\varphi^{f}\right)\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\divg\Delta\mathbf{u}\grad\varphi^{f}+\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\varphi^{s}\left(\divg\Delta\mathbf{u}\right)\dot{\mathbf{v}}^{s}+\varphi^{f}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\frac{\varphi^{s}}{\varphi^{f}}\left[\left(\left(-\frac{1}{\varphi^{f}}\right)\frac{\dot{J}}{J}+\frac{\gamma}{\beta\Delta t}\right)\divg\Delta\mathbf{u}-\grad^{T}\Delta\mathbf{u}:\mathbf{L}^{s}\right]\mathbf{w}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\divg\Delta\mathbf{u}\grad\mathbf{w}+\frac{\gamma}{\beta\Delta t}\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f^{3}}}\left[\left(\left(\varphi^{s}+1\right)\frac{1}{\varphi^{f}}\grad\varphi^{f}\divg\Delta\mathbf{u}-\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)\cdot\mathbf{w}\right]\mathbf{w}dv \end{aligned} \label{eq:Wint-lin-u-u-BFSI} \end{equation}$$ $$\begin{aligned} & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\mathbf{L}^{f}\cdot\grad\Delta\mathbf{u}\cdot\mathbf{w}+\varphi^{s}\left(\divg\Delta\mathbf{u}\right)\mathbf{a}^{f}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\rho^{s}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}\Delta\mathbf{u}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\left(-2\left(\divg\Delta\mathbf{u}\right)\mathbf{k}^{-1}+\grad^{T}\Delta\mathbf{u}\cdot\mathbf{k}^{-1}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\left(\mathbf{k}^{-1}\cdot\grad\Delta\mathbf{u}+\left(\mathbf{k}^{-1}\underbar{\otimes}\mathbf{k}^{-1}\right):\mathsf{k}:\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}dv, \end{aligned}$$ for the fluid flux equations $$\begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{u}\right]= & \int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(-\frac{\varphi^{s}}{\varphi^{f}}\left[\left(-\frac{1}{\varphi^{f}}\frac{\dot{J}}{J}+\frac{\gamma}{\beta\Delta t}\right)\divg\Delta\mathbf{u}-\grad^{T}\Delta\mathbf{u}:\mathbf{L}^{s}\right]\mathbf{w}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\varphi^{s}\left(\divg\Delta\mathbf{u}\right)\dot{\mathbf{v}}^{s}+\varphi^{f}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\divg\Delta\mathbf{u}\grad\mathbf{w}+\frac{\gamma}{\beta\Delta t}\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\frac{1}{\varphi^{f^{2}}}\left[\left(\left(\varphi^{s}+1\right)\frac{1}{\varphi^{f}}\grad\varphi^{f}\left(\divg\Delta\mathbf{u}\right)-\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)\cdot\mathbf{w}\right]\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\mathbf{L}^{f}\cdot\grad\Delta\mathbf{u}\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\left(1-\frac{\varphi^{s}}{\varphi^{f}}\right)\left(\divg\Delta\mathbf{u}\right)\rho_{T}^{f}\mathbf{a}^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\left(\left(\divg\Delta\mathbf{u}\right)\mathbf{I}-\grad^{T}\Delta\mathbf{u}\right)\cdot\grad pdv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\left(2\left(\divg\Delta\mathbf{u}\right)\mathbf{k}^{-1}-\grad^{T}\Delta\mathbf{u}\cdot\mathbf{k}^{-1}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\left(\mathbf{k}^{-1}\cdot\grad\Delta\mathbf{u}+\left(\mathbf{k}^{-1}\underbar{\otimes}\mathbf{k}^{-1}\right):\mathsf{k}:\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\left(\divg\Delta\mathbf{u}\right)\mathbf{I}+\grad^{T}\Delta\mathbf{u}\right)\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\left(\divg\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\frac{1}{\varphi^{f}}\boldsymbol{\mathcal{C}}^{\tau}:\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{D}^{w}+\mathbf{M}\cdot\grad\Delta\mathbf{u}\right)\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\frac{1}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:2\frac{\varphi^{s}}{\varphi^{f}}\left(\divg\Delta\mathbf{u}\right)\mathbf{w}\otimes\grad\varphi^{f}\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{1}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\grad^{T}\Delta\mathbf{u}\cdot\grad\varphi^{f}\right)\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{1}{\varphi^{f^{3}}}\boldsymbol{\mathcal{C}}^{\tau}:\mathbf{w}\otimes\left(-\left(\divg\Delta\mathbf{u}\right)\grad\varphi^{f}+\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\alpha_{f}\boldsymbol{\tau}:\grad\delta\mathbf{w}\cdot\left(\left(1-\frac{\varphi^{s}}{\varphi^{f}}\right)\left(\divg\Delta\mathbf{u}\right)\mathbf{I}-\grad\Delta\mathbf{u}\right)dv \end{aligned} \label{eq:Wint-lin-w-u-BFSI} \end{equation}$$ $$\begin{aligned} & +\int_{\Omega^{b}}\alpha_{f}\grad\delta\mathbf{w}:\boldsymbol{\mathcal{C}}^{\tau}:\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{D}^{w}+\mathbf{M}\cdot\grad\Delta\mathbf{u}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\grad\delta\mathbf{w}:\boldsymbol{\mathcal{C}}^{\tau}:2\frac{\varphi^{s}}{\varphi^{f^{3}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{w}\otimes\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\grad\delta\mathbf{w}:\boldsymbol{\mathcal{C}}^{\tau}:\frac{1}{\varphi^{f^{2}}}\mathbf{w}\otimes\left(-\grad^{T}\Delta\mathbf{u}\cdot\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\grad\delta\mathbf{w}:\boldsymbol{\mathcal{C}}^{\tau}:\frac{1}{\varphi^{f^{2}}}\mathbf{w}\otimes\left(-\left(\divg\Delta\mathbf{u}\right)\grad\varphi^{f}+\varphi^{s}\grad\left(\divg\Delta\mathbf{u}\right)\right)dv, \end{aligned}$$ and for the fluid dilatation equations $$\begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{u}\right]= & \int_{\Omega^{b}}\alpha_{f}\grad\delta J^{f}\cdot\left(\divg\Delta\mathbf{u}\mathbf{I}-\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta J^{f}\alpha_{f}\frac{1}{J^{f}}\left(\dot{J}^{f}\divg\Delta\mathbf{u}+\grad J^{f}\cdot\left(\left(\divg\Delta\mathbf{u}\right)\mathbf{I}-\grad\Delta\mathbf{u}\right)\cdot\mathbf{w}\right)dv\\ & +\int_{\Omega^{b}}-\delta J^{f}\alpha_{f}\left(\divg\mathbf{v}^{s}+\frac{\gamma}{\beta\Delta t}\right)\divg\Delta\mathbf{u}+\grad^{T}\Delta\mathbf{u}:\mathbf{L}^{s}dv, \end{aligned} \label{eq:Wint-lin-J-u-BFSI} \end{equation}$$ where $\mathbf{M}=\frac{\gamma}{\beta\Delta t}\mathbf{I}-\mathbf{L}^{s}$ , and $\boldsymbol{\mathcal{C}}$ is the fourth-order spatial elasticity tensor associated with $\boldsymbol{\sigma}^{e}$ [23]. As done in the CFD formulation (Section 3.5.2↑) [76], we have introduced the fourth-order tensor $\boldsymbol{\mathcal{C}}^{\tau}$ representing the tangent of the viscous stress with respect to the rate of deformation, $$\begin{equation} \boldsymbol{\mathcal{C}}^{\tau}=\frac{\partial\boldsymbol{\tau}}{\partial\mathbf{D}^{f}}\,,\label{eq:visc-stress-tan-D-1} \end{equation}$$ where $\mathbf{D}^{f}$ is the fluid rate of deformation tensor (the symmetric part of $\mathbf{L}^{f}$ ). The tensors $\boldsymbol{\mathcal{C}}_{d}$ and $\boldsymbol{\mathcal{C}}$ depend on the choice of constitutive relations for the fluid and solid constituents of $\Omega^{f}$ , respectively. In addition, $\mathsf{k}$ is the spatial fourth order permeability tensor with respect to right Cauchy-Green tensor $\mathbf{C}$ . The linearized equations include the generalized-$\alpha$ parameters because the virtual work is evaluated at the intermediate time step, while the increment itself ($\Delta\mathbf{u}$ in this case) is at the current time step [77].

Following the same procedure, the linearizations of $\delta W_{int}$ along an increment $\Delta\mathbf{w}$ is given by $$\begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{w}\right]= & \int_{\Omega^{b}}-\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}^{\tau}:\left(\grad\Delta\mathbf{w}-\frac{1}{\varphi^{f}}\Delta\mathbf{w}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\left(\frac{\alpha_{m}}{\gamma\alpha_{f}\Delta t}-\frac{\varphi^{s}}{\varphi^{f}}\frac{\dot{J}}{J}-\frac{1}{\varphi^{f^{2}}}\left(\grad\varphi^{f}\cdot\mathbf{w}\right)\right)\mathbf{I}+\mathbf{L}^{f}\right)\cdot\Delta\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f^{2}}}\grad\Delta\mathbf{w}\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\frac{\left(\varphi^{s}\right)^{2}}{\left(\varphi^{f}\right)^{2}}\boldsymbol{\mathcal{C}}^{\tau}:\left(\grad\Delta\mathbf{w}-\frac{1}{\varphi^{f}}\Delta\mathbf{w}\otimes\grad\varphi^{f}\right)\cdot\grad\left(\frac{\varphi^{f}}{\varphi^{s}}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot-\alpha_{f}\mathbf{k}^{-1}\cdot\Delta\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(\left(\frac{\alpha_{m}}{\gamma\alpha_{f}\Delta t}-\frac{\varphi^{s}}{\varphi^{f}}\frac{\dot{J}}{J}-\frac{1}{\varphi^{f^{2}}}\left(\grad\varphi^{f}\cdot\mathbf{w}\right)\right)\mathbf{I}+\mathbf{L}^{f}\right)\cdot\Delta\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f^{2}}}\grad\Delta\mathbf{w}\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\mathbf{k}^{-1}\cdot\Delta\mathbf{w}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\frac{1}{\varphi^{f^{2}}}\boldsymbol{\mathcal{C}}^{\tau}:\left(\frac{1}{\varphi^{f}}\Delta\mathbf{w}\otimes\grad\varphi^{f}-\grad\Delta\mathbf{w}\right)\cdot\grad\varphi^{f}dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{1}{\varphi^{f}}\grad\delta\mathbf{w}:\boldsymbol{\mathcal{C}}^{\tau}:\left(\grad\Delta\mathbf{w}-\frac{1}{\varphi^{f}}\Delta\mathbf{w}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\grad\delta J^{f}\cdot\Delta\mathbf{w}+\delta J^{f}\alpha_{f}\frac{1}{J^{f}}\grad J^{f}\cdot\Delta\mathbf{w}dv, \end{aligned} \label{eq:Wint-lin-w-BFSI} \end{equation}$$ whereas the linearizations of $\delta W_{int}$ along an increment $\Delta J^{f}$ is $$\begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta J^{f}\right]= & \int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\boldsymbol{\tau}_{J}^{\prime}:\grad\delta\mathbf{v}\Delta J^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\varphi^{s}\left(\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}_{J}^{\prime}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}+\frac{\rho_{T}^{f}}{J^{f}}\mathbf{a}^{f}\right)\Delta J^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\left(p^{\prime\prime}\grad J^{f}\Delta J^{f}+p^{\prime}\grad\Delta J^{f}\right)dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\left(\frac{\rho_{T}^{f}}{J^{f}}\mathbf{a}^{f}\Delta J^{f}+\frac{1}{\varphi^{f}}\boldsymbol{\tau}_{J}^{\prime}\Delta J^{f}\cdot\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\boldsymbol{\tau}_{J}^{\prime}:\grad\delta\mathbf{w}\Delta J^{f}dv\\ & +\int_{\Omega^{b}}\delta J^{f}\alpha_{f}\frac{1}{J^{f}}\left(\varphi^{f}\frac{\alpha_{m}}{\gamma\alpha_{f}\Delta t}-\frac{1}{J^{f}}\left(\varphi^{f}\dot{J}^{f}+\grad J^{f}\cdot\mathbf{w}\right)\right)\Delta J^{f}dv\\ & +\int_{\Omega^{b}}\delta J^{f}\alpha_{f}\frac{1}{J^{f}}\grad\Delta J^{f}\cdot\mathbf{w}dv. \end{aligned} \label{eq:Wint-lin-J-BFSI} \end{equation}$$ Here, $p^{\prime}$ and $p^{\prime\prime}$ respectively represent the first and second derivatives of $p\left(J^{f}\right)$ . We have also defined $\boldsymbol{\tau}_{J}^{\prime}$ as the tangent of the viscous stress $\boldsymbol{\tau}$ with respect to $J^{f}$ , $$\begin{equation} \boldsymbol{\tau}_{J}^{\prime}=\frac{\partial\boldsymbol{\tau}}{\partial J^{f}}\,.\label{eq:visc-stress-tan-J-1} \end{equation}$$

For the external work, when $\mathbf{t}^{\sigma}$ , $\mathbf{t}^{\tau}$ , all $\mathbf{b}$ and $w_{n}$ are prescribed, the linearizations simplify to $$\begin{equation} D\delta W_{ext}\left[\Delta\mathbf{u}\right]=\int_{\Omega^{b}}\delta\mathbf{w}\cdot\alpha_{f}\divg\Delta\mathbf{u}\rho_{T}^{f}\mathbf{b}^{f}dv\,,\label{eq:Wext-lin-u-BFSI} \end{equation}$$ $$\begin{equation} D\delta W_{ext}\left[\Delta\mathbf{w}\right]=0\,,\label{eq:Wext-lin-w-BFSI} \end{equation}$$

$$\begin{equation} \begin{aligned}D\delta W_{ext}\left[\Delta J^{f}\right]= & \int_{\Omega^{b}}\delta\mathbf{w}\cdot-\alpha_{f}\frac{\rho_{T}^{f}}{J^{f}}\mathbf{b}^{f}\Delta J^{f}dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{J^{f}}\mathbf{b}^{f}\Delta J^{f}dv\,. \end{aligned} \label{eq:Wext-lin-J-BFSI} \end{equation}$$

As discussed in Section ↓, we define the fluid dilatation $e^{f}=J^{f}-1$ as an alternative essential variable, since initial and boundary conditions $e^{f}=0$ are more convenient to handle in a numerical scheme than $J^{f}=1$ . It follows that $\grad J^{f}=\grad e^{f}$ and $\dot{J}^{f}=\dot{e}^{f}$ . Therefore the changes to the above equations are minimal, simply requiring the substitution $J^{f}=1+e^{f}$ and $\Delta J^{f}=\Delta e^{f}$ .