Virtual Work and Weak Form
The virtual work statement is used to enforce the three governing equations needed to solve for the nodal DOFs \(\mathbf{u}\), \(\mathbf{w}\) and \(e^{f}\), namely the mixture mass balance (2.15-6), the fluid momentum balance (2.15-7), and the solid momentum balance (2.15-8). We may rewrite the momentum balance equations to facilitate the enforcement of natural traction boundary conditions given in (2.15-16) and (2.15-17). Using \(\boldsymbol{\tau}^{f}=\varphi^{f}\boldsymbol{\tau}\), these become
\[ \begin{equation} \begin{aligned}\varphi^{s}\left(-\rho_{T}^{f}\mathbf{a}^{f}+\rho_{T}^{s}\mathbf{a}^{s}\right)= & -\frac{\varphi^{s^{2}}}{\varphi^{f}}\mathbf{\boldsymbol{\tau}}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}+\divg\left(-\varphi^{s}\boldsymbol{\tau}+\boldsymbol{\sigma}^{e}\right)\\ & +\varphi^{s}\left(-\rho_{T}^{f}\mathbf{b}^{f}+\rho_{T}^{s}\mathbf{b}^{s}\right)+\mathbf{k}^{-1}\cdot\mathbf{w}\,,\\ \rho_{T}^{f}\mathbf{a}^{f}= & -\grad p+\frac{1}{\phi^{f}}\mathbf{\boldsymbol{\tau}}\cdot\grad\varphi^{f}+\divg\boldsymbol{\tau}\\ & +\rho_{T}^{f}\mathbf{b}^{f}-\mathbf{k}^{-1}\cdot\mathbf{w}\,. \end{aligned} \label{eq:bfsi-gov-eqn-jump} \end{equation} \]
The virtual work statement for a Galerkin finite element formulation is \(\delta W=0\), where
\[ \begin{equation} \begin{aligned}\delta W & =\int_{\Omega^{b}}\delta\mathbf{v}^{s}\cdot\left[\divg\left(-\varphi^{s}\boldsymbol{\tau}+\boldsymbol{\sigma}^{e}\right)-\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}+\mathbf{k}^{-1}\cdot\mathbf{w}\right]\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}^{s}\cdot\varphi^{s}\left(-\rho_{T}^{f}\left(\mathbf{b}^{f}-\mathbf{a}^{f}\right)+\rho_{T}^{s}\left(\mathbf{b}^{s}-\mathbf{a}^{s}\right)\right)\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\left[\divg\boldsymbol{\tau}-\grad p+\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\phi^{f}-\mathbf{k}^{-1}\cdot\mathbf{w}\right]\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\rho_{T}^{f}\left(\mathbf{b}^{f}-\mathbf{a}^{f}\right)\,dv\\ & +\int_{\Omega^{b}}\delta J^{f}\left[\divg\mathbf{w}+\frac{\dot{J}^{s}}{J^{s}}-\frac{1}{J^{f}}\left(\varphi^{f}\dot{J}^{f}+\grad J^{f}\cdot\mathbf{w}\right)\right]\,dv\,. \end{aligned} \label{eq:bfsi-virtual-work-strong} \end{equation} \]
These integrals are evaluated in the current configuration of \(\Omega^{b}\). Here, \(\delta\mathbf{v}^{s}\) is the virtual solid velocity, \(\delta\mathbf{w}\) is the virtual relative fluid volumetric flux, and \(\delta J^{f}\) is the virtual fluid energy density. Integrating by parts and using the divergence theorem, the weak form of this statement may be written as \(\delta W=\delta W_{ext}-\delta W_{int}\) where the internal virtual work is
\[ \begin{equation} \begin{aligned}\delta W_{int} & =\int_{\Omega^{b}}\left[\left(-\varphi^{s}\boldsymbol{\tau}+\boldsymbol{\sigma}^{e}\right):\grad\delta\mathbf{v}^{s}-\delta\mathbf{v}^{s}\cdot\mathbf{k}^{-1}\cdot\mathbf{w}\right]\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{v}^{s}\cdot\varphi^{s}\left(\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}-\rho_{T}^{f}\mathbf{a}^{f}+\rho_{T}^{s}\mathbf{a}^{s}\right)\,dv\\ & +\int_{\Omega^{b}}\boldsymbol{\tau}:\grad\delta\mathbf{w}\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\left(\grad p+\mathbf{k}^{-1}\cdot\mathbf{w}\right)\,dv\\ & +\int_{\Omega^{b}}\delta\mathbf{w}\cdot\left(\rho_{T}^{f}\mathbf{a}^{f}-\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\varphi^{f}\right)\,dv\\ & +\int_{\Omega^{b}}\mathbf{w}\cdot\grad\delta J^{f}\,dv\,,\\ & +\int_{\Omega^{b}}\delta J^{f}\left[\frac{1}{J^{f}}\left(\varphi^{f}\dot{J}^{f}+\grad J^{f}\cdot\mathbf{w}\right)-\frac{\dot{J}^{s}}{J^{s}}\right]\,dv \end{aligned} \label{eq:bfsi-int-virtual-work} \end{equation} \]
and the external part is
\[ \begin{equation} \begin{aligned}\delta W_{ext} & =\int_{\partial\Omega^{b}}\delta\mathbf{v}^{s}\cdot\mathbf{t}^{\sigma}\,da\\ & +\int_{\Omega^{b}}\delta\mathbf{v}^{s}\cdot\varphi^{s}\left(-\rho_{T}^{f}\mathbf{b}^{f}+\rho_{T}^{s}\mathbf{b}^{s}\right)\,dv\\ & +\int_{\partial\Omega^{b}}\delta\mathbf{w}\cdot\mathbf{t}^{\tau}\,da+\int_{\Omega^{f}}\delta\mathbf{w}\cdot\rho_{T}^{f}\mathbf{b}^{f}\,dv\\ & +\int_{\partial\Omega^{b}}\delta J^{f}w_{n}\,da\,. \end{aligned} \label{eq:bfsi-ext-virtual-work} \end{equation} \]
where
\[ \begin{equation} \mathbf{t}^{\sigma}=-\varphi^{s}\mathbf{t}^{\tau}+\mathbf{t}^{e}\,.\label{eq:bfsi-sigma-traction} \end{equation} \]
Here, \(\mathbf{t}^{e}=\boldsymbol{\sigma}^{e}\cdot\mathbf{n}\) is the elastic traction, \(\mathbf{t}^{\tau}=\boldsymbol{\tau}\cdot\mathbf{n}\) is the true fluid viscous traction, and \(w_{n}=\mathbf{w}_{n}\cdot\mathbf{n}\) is the normal component of the relative fluid flux on the boundary \(\partial\Omega^{b}\), whose outward unit normal is \(\mathbf{n}\). The traction \(\mathbf{t}^{\sigma}\) emerges from the jump condition in (2.15-17). The integrands of the surface integrals represent the natural boundary conditions for this formulation. If boundary conditions are not set explicitly on \(\partial\Omega^{b}\), the natural boundary conditions are \(\mathbf{t}^{\sigma}=\mathbf{0}\), \(\mathbf{t}^{\tau}=\mathbf{0}\), and \(w_{n}=0\). These natural boundary conditions are consistent with the jump conditions presented above. Essential boundary conditions are prescribed on the solid displacement \(\mathbf{u}\), relative volumetric fluid flux \(\mathbf{w}\), and fluid dilatation \(e^{f}\), which are also consistent with the above jump conditions. In particular, an essential no-slip boundary condition may be prescribed on \(\Gamma^{bs}\) by setting \(\mathbf{w}=\mathbf{0}\). A symmetry plane may be prescribed with the essential boundary condition \(u_{n}\equiv\mathbf{u}\cdot\mathbf{n}=0\) and the natural boundary conditions \(\mathbf{t}^{\tau}=\mathbf{0}\) and \(w_{n}=0\).
In this formulation, the mixture traction is defined as \(\mathbf{t}=-p\mathbf{n}+\mathbf{t}^{e}+\varphi^{f}\mathbf{t}^{\tau}\), which may also be written as \(\mathbf{t}=-p\mathbf{n}+\mathbf{t}^{\sigma}+\mathbf{t}^{\tau}\). Because of the way we chose to split the internal and external virtual work in \eqref{eq:bfsi-int-virtual-work}-\eqref{eq:bfsi-ext-virtual-work}, \(\mathbf{t}\) is not a natural boundary condition in this formulation. In this expression for \(\mathbf{t}\), \(\mathbf{t}^{\sigma}\), and \(\mathbf{t}^{\tau}\) may be prescribed as natural boundary conditions, whereas \(p\) may be prescribed as an essential boundary condition on \(e^{f}\), using (2.15-9). However, there are two general scenarios where \(\mathbf{t}\) needs to be prescribed on a region of the boundary \(\partial\Omega^{b}\) with incomplete prior knowledge of \(p\), \(\mathbf{t}^{\sigma}\), or \(\mathbf{t}^{\tau}\): (1) When a BFSI boundary \(\Gamma^{b}\) represents a free surface (such as the fluid surface in channel flow), the mixture traction boundary condition requires that \(\mathbf{t}=\mathbf{0}\), in which case it is necessary to explicitly enforce \(-p\mathbf{n}+\mathbf{t}^{\sigma}+\mathbf{t}^{\tau}=\mathbf{0}\) as a constraint equation on that boundary, to impart the free surface \(\Gamma^{b}\) its natural shape. (2) At a biphasic-solid interface \(\Gamma^{bs}\), \(\mathbf{t}\) must balance the traction acting on the adjoining solid domain. Since \(\mathbf{u}\) is continuous across \(\Gamma^{bs}\) due to shared nodes, the solid natural boundary condition \(\mathbf{t}^{\sigma}\) of \(\mathbf{t}\) is already accounted for by the deformation, so that it is only necessary to prescribe the portion \(-p\mathbf{n}+\mathbf{t}^{\tau}\) of \(\mathbf{t}\) on the solid domain, thus \(-p\mathbf{n}+\mathbf{t}^{\tau}+\mathbf{t}^{s}=\mathbf{0}\) where \(\mathbf{t}^{s}\) is the (equal and opposite) traction acting on the solid domain. In both cases, the form of this traction boundary condition is the same, with \(\mathbf{t}^{e}\) and \(\mathbf{t}^{s}\) representing the tractions acting on \(\Gamma^{b}\) and \(\Gamma^{bs}\), respectively.
For both of these cases, the resulting virtual work on the free surface \(\Gamma^{b}\) or the interface \(\Gamma^{bs}\) takes the form
\[ \begin{equation} \delta F=-\int_{\Gamma}\delta\mathbf{v}^{s}\cdot\left(-p\mathbf{n}+\mathbf{t}^{\tau}\right)\,da\,,\label{eq:bfsi-traction-virtual-work} \end{equation} \]
where the elemental area \(da\) on \(\Gamma\) may be evaluated from the covariant basis vectors \(\mathbf{g}_{\alpha}\) (\(\alpha=1,2\)),
\[ \begin{equation} da=\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\,d\eta^{1}d\eta^{2}\,,\label{eq:bfsi-diff-area} \end{equation} \]
where
\[ \begin{equation} \mathbf{g}_{\alpha}=\frac{\partial\mathbf{x}\left(\eta^{1},\eta^{2}\right)}{\partial\eta^{\alpha}}\,,\label{eq:bfsi-covar-vectors} \end{equation} \]
and \(\mathbf{x}\left(\eta^{1},\eta^{2}\right)\) is the parametric representation of \(\Gamma\), defined on the solid constituent. The outward normal \(\mathbf{n}\) to \(\Omega^{b}\) on \(\Gamma\) is evaluated from
\[ \begin{equation} \mathbf{n}=\frac{\mathbf{g}_{1}\times\mathbf{g}_{2}}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,.\label{eq:bfsi-norm-vector} \end{equation} \]
As a result, the virtual work can be rewritten as
\[ \begin{equation} \delta F=-\int_{\Gamma}\delta\mathbf{v}^{s}\cdot\left(-p\mathbf{I}+\boldsymbol{\tau}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\,d\eta^{1}d\eta^{2}\,.\label{eq:bfsi-traction-final} \end{equation} \]
In FEBio this boundary condition is called BFSI traction, which the user must explicitly prescribe on free surfaces \(\Gamma^{b}\) and deformable interfaces \(\Gamma^{bs}\). The code automatically determines which of these two types of interfaces is being considered.
BFSI Linearization
Using the virtual work integral \(\delta W\) such that
\[ \begin{equation} \delta W+D\delta W\left[\Delta\mathbf{u}\right]+D\delta W\left[\Delta\mathbf{w}\right]+D\delta W\left[\Delta J^{f}\right]\approx0\,,\label{eq:Virtual-Work-Expanded-BFSI} \end{equation} \]
it may be expanded as
\[ \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. \eqref{eq:bfsi-int-virtual-work}-\eqref{eq:bfsi-ext-virtual-work}. 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\thinspace, \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 . 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\thinspace, \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\thinspace, \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\thinspace, \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}\) . As done in the CFD formulation (Section Temporal Discretization and Linearization) , 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 .
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\thinspace, \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\thinspace. \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 Mass and Momentum Balance, 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}\).
BFSI Spatial Discretization
The degrees of freedom \(\mathbf{u}\left(\mathbf{x},t\right)\), \(\mathbf{w}\left(\mathbf{x},t\right)\), \(J^{f}\left(\mathbf{x},t\right)\) are spatially interpolated over the domain \(\Omega^{f}\) using the same interpolation functions \(N_{a}\left(\mathbf{x}\right)\), with \(a=1\) to \(n\), where \(n\) is the number of nodes in an element,
\[ \begin{equation} \begin{aligned}\mathbf{u}\left(\mathbf{x},t\right) & =\sum_{a=1}^{n}N_{a}\left(\mathbf{x}\right)\mathbf{u}_{a}\,,\\ \mathbf{w}\left(\mathbf{x},t\right) & =\sum_{a=1}^{n}N_{a}\left(\mathbf{x}\right)\mathbf{w}_{a}\,,\\ J^{f}\left(\mathbf{x},t\right) & =\sum_{a=1}^{n}N_{a}\left(\mathbf{x}\right)J_{a}^{f}\,. \end{aligned} \label{eq:spatial-interpol-BFSI} \end{equation} \]
Here, \(\mathbf{u}_{a}\), \(\mathbf{w}_{a}\), and \(J_{a}^{f}\) are the nodal values of the degrees of freedom that evolve over time. These relations may be used to evaluate \(\mathbf{L}^{s}\), \(\mathbf{L}^{w}\), \(\divg\mathbf{w}\), \(\dot{\mathbf{w}}\), \(\grad J^{f}\), \(\dot{J^{f}}\), etc. Similar interpolations are used for virtual increments \(\delta\mathbf{v}\), \(\delta\mathbf{w}\) and \(\delta J^{f}\), as well as real increments \(\Delta\mathbf{u}\), \(\Delta\mathbf{w}\) and \(\Delta J^{f}\). In practice, interpolations \(N_{a}\) are performed in the parametric space of each finite element, which is a material frame.
Now, the discretized form of \(\delta W_{int}\), using eq. \eqref{eq:bfsi-int-virtual-work}, may be written as
\[ \begin{equation} \begin{aligned}\delta W_{int} & =\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}^{u}\\ & +\sum_{a}\delta\mathbf{w}_{a}\cdot\mathbf{f}_{a}^{w}\\ & +\sum_{a}\delta J_{a}^{f}f_{a}^{J}\,, \end{aligned} \label{eq:Wint-discret-BFSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{f}_{a}^{u}= & \int_{\Omega^{b}}\left(-\varphi^{s}\boldsymbol{\tau}+\boldsymbol{\sigma}^{e}\right)\cdot\grad N_{a}dv\\ & +\int_{\Omega^{b}}N_{a}\left(\varphi^{s}\left(\frac{\varphi^{s}}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\frac{\varphi^{f}}{\varphi^{s}}-\rho_{T}^{f}\mathbf{a}^{f}+\rho_{T}^{s}\dot{\mathbf{v}}^{s}\right)-\mathbf{k}^{-1}\cdot\mathbf{w}\right)dv\,,\\ \mathbf{f}_{a}^{w}= & \int_{\Omega^{b}}\boldsymbol{\tau}\cdot\grad N_{a}+N_{a}\left(\rho_{T}^{f}\mathbf{a}^{f}+\grad p-\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\grad\varphi^{f}+\mathbf{k}^{-1}\cdot\mathbf{w}\right)dv\,,\\ f_{a}^{J}= & \int_{\Omega^{b}}\grad N_{a}\cdot\mathbf{w}+N_{a}\left(\frac{\varphi^{f}}{J^{f}}\frac{D^{f}J^{f}}{Dt}-\frac{\dot{J}^{s}}{J^{s}}\right)dv\,. \end{aligned} \label{eq:Wint-discret-parts-BFSI} \end{equation} \]
The discretized form of \(D\delta W_{int}\left[\Delta\mathbf{u}\right]\) becomes
\[ \begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{u}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{uu}\cdot\Delta\mathbf{u}_{b}\\ & +\sum_{a}\delta\mathbf{w}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{wu}\cdot\Delta\mathbf{u}_{b}\\ & +\sum_{a}\delta J_{a}^{f}\sum_{b}\mathbf{k}_{ab}^{Ju}\cdot\Delta\mathbf{u}_{b}\,, \end{aligned} \label{eq:Wint-lin-u-discret-FSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{uu}= & \int_{\Omega^{b}}\alpha_{f}\left(\left(\boldsymbol{\sigma}^{e}:\grad N_{b}\otimes\grad N_{a}\right)\mathbf{I}+\grad N_{a}\cdot\boldsymbol{\mathcal{C}}\cdot\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\varphi^{s}\boldsymbol{\tau}\cdot\left(\frac{\varphi^{s}}{\varphi^{f}}\grad N_{a}\otimes\grad N_{b}+\grad N_{b}\otimes\grad N_{a}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\grad N_{b}\cdot\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{D}^{w}+\mathbf{M}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{\varphi^{s}}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\varphi^{f}\otimes\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{\varphi^{s}}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(\grad N_{b}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{\varphi^{s}}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(-\grad\varphi^{f}\otimes\grad N_{b}+\varphi^{s}\grad\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\tau}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\frac{\varphi^{f}}{\varphi^{s}}\otimes\grad N_{b}+\grad N_{b}\otimes\grad\frac{\varphi^{f}}{\varphi^{s}}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f}}\boldsymbol{\tau}\cdot\frac{1}{\varphi^{s}}\grad\grad N_{b}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f}}\grad\frac{\varphi^{f}}{\varphi^{s}}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\grad N_{b}\cdot\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{D}^{w}+\mathbf{M}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\grad\frac{\varphi^{f}}{\varphi^{s}}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\varphi^{f}\otimes\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\grad\frac{\varphi^{f}}{\varphi^{s}}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(\grad N_{b}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\grad\frac{\varphi^{f}}{\varphi^{s}}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(-\grad\varphi^{f}\otimes\grad N_{b}+\varphi^{s}\grad\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\rho_{T}^{f}N_{a}\left(\varphi^{s}\dot{\mathbf{v}}^{s}\otimes\grad N_{b}+\varphi^{f}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}N_{b}\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\rho_{T}^{f}N_{a}\frac{\varphi^{s^{2}}}{\varphi^{f^{2}}}\left[\left(\left(-\frac{1}{\varphi^{f}}\right)\frac{\dot{J}}{J}+\frac{\gamma}{\beta\Delta t}\right)\mathbf{w}\otimes\grad N_{b}-\mathbf{w}\otimes\left(\mathbf{L}^{s^{T}}\cdot\grad N_{b}\right)\right]dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{L}^{w}\cdot\mathbf{w}\otimes\grad N_{b}+\frac{\gamma}{\beta\Delta t}\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\rho_{T}^{f}\frac{\varphi^{s}}{\varphi^{f^{4}}}\left(\varphi^{s}+1\right)\left(\grad\varphi^{f}\cdot\mathbf{w}\right)\mathbf{w}\otimes\grad N_{b}dv \end{aligned} \label{eq:Wint-lin-uu-discret-BFSI} \end{equation} \]
\[ \begin{aligned} & +\int_{\Omega^{b}}\alpha_{f}N_{a}\rho_{T}^{f}\frac{\varphi^{s^{2}}}{\varphi^{f^{3}}}\mathbf{w}\otimes\grad\grad^{T}N_{b}\cdot\mathbf{w}dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\rho_{T}^{f}N_{a}\left(\mathbf{L}^{f}\left(\grad N_{b}\cdot\mathbf{w}\right)+\varphi^{s}\mathbf{a}^{f}\otimes\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}\rho^{s}N_{b}\mathbf{I}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(-2\left(\mathbf{k}^{-1}\cdot\mathbf{w}\right)\otimes\grad N_{b}+\grad N_{b}\otimes\left(\mathbf{k}^{-1}\cdot\mathbf{w}\right)\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{k}^{-1}+\left(\mathbf{k}^{-1}\underbar{\otimes}\mathbf{k}^{-1}\right):\mathsf{k}:\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}\right)dv\,, \end{aligned} \]
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{wu}= & \int_{\Omega^{b}}\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}N_{a}\left(\varphi^{s}\dot{\mathbf{v}}^{s}\otimes\grad N_{b}+\varphi^{f}\frac{\alpha_{m}}{\alpha_{f}\beta\Delta t^{2}}N_{b}\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}N_{a}\frac{\varphi^{s}}{\varphi^{f}}\left(-\frac{1}{\varphi^{f}}\frac{\dot{J}}{J}+\frac{\gamma}{\beta\Delta t}\right)\mathbf{w}\otimes\grad N_{b}dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}N_{a}\frac{\varphi^{s}}{\varphi^{f}}\mathbf{w}\otimes\left(\mathbf{L}^{s^{T}}\cdot\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\rho_{T}^{f}}{\varphi^{f}}\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(\mathbf{L}^{w}\cdot\mathbf{w}\right)\otimes\grad N_{b}+\frac{\gamma}{\beta\Delta t}\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\rho_{T}^{f}}{\varphi^{f^{4}}}\left(\varphi^{s}+1\right)\left(\grad\varphi^{f}\cdot\mathbf{w}\right)\mathbf{w}\otimes\grad N_{b}dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{\rho_{T}^{f}}{\varphi^{f^{3}}}\varphi^{s}\mathbf{w}\otimes\left(\grad\grad^{T}N_{b}\cdot\mathbf{w}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{\rho_{T}^{f}}{\varphi^{f}}N_{a}\mathbf{L}^{f}\left(\grad N_{b}\cdot\mathbf{w}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(1-\frac{\varphi^{s}}{\varphi^{f}}\right)\rho_{T}^{f}\mathbf{a}^{f}\otimes\grad N_{b}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(\grad p\otimes\grad N_{b}-\grad N_{b}\otimes\grad p\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(2\left(\mathbf{k}^{-1}\cdot\mathbf{w}\right)\otimes\grad N_{b}-\grad N_{b}\otimes\left(\mathbf{k}^{-1}\cdot\mathbf{w}\right)\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\left(\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{k}^{-1}+\left(\mathbf{k}^{-1}\underbar{\otimes}\mathbf{k}^{-1}\right):\mathsf{k}:\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\varphi^{f}\otimes\grad N_{b}+\grad N_{b}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{1}{\varphi^{f}}\boldsymbol{\tau}\cdot\varphi^{s}\grad\grad N_{b}dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{1}{\varphi^{f}}\grad\varphi^{f}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\grad N_{b}\cdot\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{D}^{w}+\mathbf{M}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{1}{\varphi^{f^{3}}}\grad\varphi^{f}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\varphi^{f}\otimes\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\frac{1}{\varphi^{f^{3}}}\grad\varphi^{f}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(\grad N_{b}\otimes\grad\varphi^{f}\right)dv \end{aligned} \label{eq:Wint-lin-wu-discret-BFSI} \end{equation} \]
\[ \begin{aligned} & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{1}{\varphi^{f^{3}}}\grad\varphi^{f}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(-\grad\varphi^{f}\otimes\grad N_{b}+\varphi^{s}\grad\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\boldsymbol{\tau}\cdot\left(\left(1-\frac{\varphi^{s}}{\varphi^{f}}\right)\grad N_{a}\otimes\grad N_{b}-\grad N_{b}\otimes\grad N_{a}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\grad N_{b}\cdot\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{D}^{w}+\mathbf{M}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{1}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f}}\grad\varphi^{f}\otimes\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{1}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(\grad N_{b}\otimes\grad\varphi^{f}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\frac{1}{\varphi^{f^{2}}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(-\grad\varphi^{f}\otimes\grad N_{b}+\varphi^{s}\grad\grad N_{b}\right)dv\,, \end{aligned} \]
\[ \begin{equation} \begin{aligned}\mathbf{k}_{ab}^{Ju}= & \int_{\Omega^{b}}\alpha_{f}\left(\grad N_{b}\otimes\mathbf{w}-\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}\right)\cdot\grad N_{a}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{1}{J^{f}}\left(\dot{J}^{f}\grad N_{b}+\left(\grad N_{b}\otimes\mathbf{w}-\grad N_{b}\cdot\mathbf{w}\mathbf{I}\right)\cdot\grad J^{f}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}\left(\frac{\dot{J}}{J}+\frac{\gamma}{\beta\Delta t}\right)\grad N_{b}+\mathbf{L}^{s^{T}}\cdot\grad N_{b}dv\,, \end{aligned} \label{eq:Wint-lin-Ju-discret-BFSI} \end{equation} \]
whereas that of \(D\delta W_{int}\left[\Delta\mathbf{w}\right]\) becomes
\[ \begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta\mathbf{w}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{uw}\cdot\Delta\mathbf{w}_{b}\\ & +\sum_{a}\delta\mathbf{w}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{ww}\cdot\Delta\mathbf{w}_{b}\\ & +\sum_{a}\delta J_{a}^{f}\sum_{b}\mathbf{k}_{ab}^{Jw}\cdot\Delta\mathbf{w}_{b}\,, \end{aligned} \label{eq:Wint-lin-w-discret-BFSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{uw}= & \int_{\Omega^{b}}\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\left(\frac{1}{\varphi^{f}}N_{b}\grad\varphi^{f}-\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f}}N_{a}\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)N_{b}dv\\ & +\int_{\Omega^{b}}-\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{\varphi^{f^{2}}}N_{a}\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\left(\varphi^{s}\right)^{2}}{\left(\varphi^{f}\right)^{2}}\grad\frac{\varphi^{f}}{\varphi^{s}}\cdot\boldsymbol{\mathcal{C}}_{d}\cdot\left(-\frac{1}{\varphi^{f}}N_{b}\grad\varphi^{f}+\grad N_{b}\right)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{a}N_{b}\mathbf{k}^{-1}dv\,,\\ \mathbf{K}_{ab}^{ww}= & \int_{\Omega^{b}}\alpha_{f}N_{a}\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)N_{b}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\frac{\rho_{T}^{f}}{\varphi^{f^{2}}}\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}dv\\ & +\int_{\Omega^{b}}\alpha_{f}N_{a}\left(\frac{1}{\varphi^{f^{2}}}\grad\varphi^{f}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\left(\frac{1}{\varphi^{f}}N_{b}\grad\varphi^{f}-\grad N_{b}\right)+N_{b}\mathbf{k}^{-1}\right)dv\\ & +\int_{\Omega^{b}}\alpha_{f}\frac{1}{\varphi^{f}}\grad N_{a}\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\left(-\frac{1}{\varphi^{f}}N_{b}\grad\varphi^{f}+\grad N_{b}\right)dv\,,\\ \mathbf{k}_{ab}^{Jw}= & \int_{\Omega^{b}}\alpha_{f}N_{b}\left(\grad N_{a}+N_{a}\frac{1}{J^{f}}\grad J^{f}\right)dv\,, \end{aligned} \label{eq:Wint-lin-w-discret-parts-BFSI} \end{equation} \]
and finally, for \(D\delta W_{int}\left[\Delta J^{f}\right]\) the equations become
\[ \begin{equation} \begin{aligned}D\delta W_{int}\left[\Delta J^{f}\right]= & \sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{uJ}\Delta J_{b}^{f}\\ & +\sum_{a}\delta\mathbf{w}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{wJ}\Delta J_{b}^{f}\\ & +\sum_{a}\delta J_{a}^{f}\sum_{b}k_{ab}^{JJ}\Delta J_{b}^{f}\,, \end{aligned} \label{eq:Wint-lin-J-discret-BFSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{k}_{ab}^{uJ}= & \int_{\Omega^{b}}\alpha_{f}\varphi^{s}N_{a}N_{b}\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)dv\\ & +\int_{\Omega^{b}}-\alpha_{f}N_{b}\varphi^{s}\boldsymbol{\tau}_{J}^{\prime}\cdot\grad N_{a}dv\,,\\ \mathbf{k}_{ab}^{wJ}= & \int_{\Omega^{b}}\alpha_{f}N_{a}\left(p^{\prime}\grad N_{b}+N_{b}\left(p^{\prime\prime}\grad J^{f}-\frac{\rho_{T}^{f}}{J^{f}}\mathbf{a}^{f}-\frac{1}{\varphi^{f}}\boldsymbol{\tau}_{J}^{\prime}\cdot\grad\varphi^{f}\right)\right)dv\Delta J_{b}^{f}\\ & +\int_{\Omega^{b}}\alpha_{f}N_{b}\boldsymbol{\tau}_{J}^{\prime}\cdot\grad N_{a}dv\,,\\ k_{ab}^{JJ}= & \int_{\Omega^{b}}\alpha_{f}N_{a}\frac{1}{J^{f}}\left(\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)N_{b}+\grad N_{b}\cdot\mathbf{w}\right)dv\,, \end{aligned} \label{eq:Wint-lin-J-discret-parts-BFSI} \end{equation} \]
for external virtual work in eq. \eqref{eq:bfsi-ext-virtual-work}, the discretized equations are
\[ \begin{equation} \begin{aligned}\delta W_{ext} & =\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}^{u}\\ & +\sum_{a}\delta\mathbf{w}_{a}\cdot\mathbf{f}_{a}^{w}\\ & +\sum_{a}\delta J_{a}^{f}f_{a}^{J}\,, \end{aligned} \label{eq:Wext-discret-BFSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{f}_{a,ext}^{u} & =\int_{\partial\Omega^{b}}N_{a}\mathbf{t}^{\sigma}da+\int_{\Omega^{b}}N_{a}\varphi^{s}\left(-\rho_{T}^{f}\mathbf{b}^{f}+\rho_{T}^{s}\mathbf{b}^{s}\right)dv\,,\\ \mathbf{f}_{a,ext}^{w} & =\int_{\partial\Omega^{b}}N_{a}\mathbf{t}^{\tau}da+\int_{\Omega^{b}}\rho_{T}^{f}N_{a}\mathbf{b}^{f}dv\,,\\ f_{a,ext}^{J} & =\int_{\partial\Omega^{b}}N_{a}w_{n}da\,, \end{aligned} \label{eq:Wext-discret-parts-BFSI} \end{equation} \]
and the discretized forms of the linearized external virtual work are
\[ \begin{equation} \begin{aligned}D\delta W_{ext}\left[\Delta\mathbf{u}\right] & =\sum_{a}\delta\mathbf{w}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{wu}\cdot\Delta\mathbf{u}_{b}\,,\\ D\delta W_{ext}\left[\Delta J^{f}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{uJ}\Delta J_{b}^{f}+\sum_{a}\delta\mathbf{w}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{wJ}\Delta J_{b}^{f}\,, \end{aligned} \label{eq:Wext-lin-discret-BFSI} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab,ext}^{wu}= & \int_{\Omega^{b}}\alpha_{f}N_{a}\rho_{T}^{f}\mathbf{b}^{f}\otimes\grad N_{b}dv\,,\\ \mathbf{k}_{ab,ext}^{uJ}= & -\int_{\Omega^{b}}\alpha_{f}\varphi^{s}\frac{\rho_{T}^{f}}{J^{f}}N_{a}N_{b}\mathbf{b}^{f}dv\,,\\ \mathbf{k}_{ab,ext}^{wJ}= & -\int_{\Omega^{b}}\alpha_{f}\frac{\rho_{T}^{f}}{J^{f}}N_{a}N_{b}\mathbf{b}^{f}dv\,. \end{aligned} \label{eq:Wext-lin-discret-parts-BFSI} \end{equation} \]
Combining these results, from the linearized virtual work equation of (3.5-5) we can represent the system of equations in a compact matrix form as
\[ \begin{equation} \left[\begin{array}{ccc} \mathbf{K}_{ab}^{uu} & \mathbf{K}_{ab}^{uw} & \mathbf{k}_{ab}^{uJ}-\mathbf{k}_{ab,ext}^{uJ}\\ \mathbf{K}_{ab}^{wu}-\mathbf{K}_{ab,ext}^{wu} & \mathbf{K}_{ab}^{ww} & \mathbf{k}_{ab}^{wJ}-\mathbf{k}_{ab,ext}^{wJ}\\ \mathbf{k}_{ab}^{Ju} & \mathbf{k}_{ab}^{Jw} & k_{ab}^{JJ} \end{array}\right]\left[\begin{array}{c} \Delta\mathbf{u}\\ \Delta\mathbf{w}\\ \Delta J^{f} \end{array}\right]=\left[\begin{array}{c} \mathbf{f}_{a,ext}^{u}-\mathbf{f}_{a}^{u}\\ \mathbf{f}_{a,ext}^{w}-\mathbf{f}_{a}^{w}\\ f_{a,ext}^{J}-f_{a}^{J} \end{array}\right]\,.\label{eq:Virtual-work-matrix-BFSI} \end{equation} \]
BFSI Traction Interface
The virtual work \(\delta F\) on a biphasic-FSI interface was given in eq. \eqref{eq:bfsi-traction-final}. The linearizations of \(\delta F\) are given by
\[ \begin{equation} \begin{aligned}D\delta F\left[\Delta\mathbf{u}\right]= & \int_{\partial\Omega^{b}}\delta\mathbf{v}\cdot\alpha_{f}\frac{\varphi^{s}}{\varphi^{f}}\left(\divg\Delta\mathbf{u}\right)\mathbf{T}_{d}\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\boldsymbol{\mathcal{C}}^{\tau}:\alpha_{f}\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\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\boldsymbol{\mathcal{C}}^{\tau}:\alpha_{f}\left(2\frac{\varphi^{s}}{\varphi^{f^{3}}}\left(\divg\Delta\mathbf{u}\right)\mathbf{w}\otimes\grad\varphi^{f}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\boldsymbol{\mathcal{C}}^{\tau}:\alpha_{f}\left(-\frac{1}{\varphi^{f^{2}}}\mathbf{w}\otimes\left(-\grad^{T}\Delta\mathbf{u}\cdot\grad\varphi^{f}\right)\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}\delta\mathbf{v}\cdot\boldsymbol{\mathcal{C}}^{\tau}:\alpha_{f}\left(\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{w}\otimes\left(-\frac{1}{J}\left(\divg\Delta\mathbf{u}\right)\grad J+\grad\left(\divg\Delta\mathbf{u}\right)\right)\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\alpha_{f}\boldsymbol{\tau}\cdot\left(\frac{\partial\Delta\mathbf{u}}{\partial\eta^{1}}\times\mathbf{g}_{2}+\mathbf{g}_{1}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{2}}\right)\thinspace d\eta^{1}d\eta^{2}\,,\\ D\delta F\left[\Delta\mathbf{w}\right] & =\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\alpha_{f}\boldsymbol{\mathcal{C}}^{\tau}:\left(\frac{1}{\varphi^{f}}\grad\Delta\mathbf{w}-\frac{1}{\varphi^{f^{2}}}\Delta\mathbf{w}\otimes\grad\varphi^{f}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\,,\\ D\delta F\left[\Delta J^{f}\right] & =\int_{\partial\Omega^{b}}-\delta\mathbf{v}\cdot\alpha_{f}\left(-p^{\prime}\mathbf{I}+\boldsymbol{\tau}_{J}^{\prime}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\Delta J^{f}\thinspace d\eta^{1}d\eta^{2}\,. \end{aligned} \label{eq:BFSI-int-lin} \end{equation} \]
The discretized forms of this virtual work is
\[ \begin{equation} \delta F=\sum_{a}\delta\mathbf{v}_{a}\cdot\mathbf{f}_{a}\,,\label{eq:BSI-int-discret} \end{equation} \]
where
\[ \begin{equation} \mathbf{f}_{a}=\int_{\partial\Omega^{b}}-N_{a}\left(-p\mathbf{I}+\boldsymbol{\tau}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\,.\label{eq:BFSI-int-discret-parts} \end{equation} \]
The discretized forms of the linearizations are
\[ \begin{equation} \begin{aligned}D\delta F\left[\Delta\mathbf{u}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{uu}\cdot\Delta\mathbf{u}_{b}\,,\\ D\delta F\left[\Delta\mathbf{w}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{K}_{ab}^{uw}\cdot\Delta\mathbf{w}_{b}\,,\\ D\delta F\left[\Delta J^{f}\right] & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\mathbf{k}_{ab}^{uJ}\Delta J_{b}^{f}\,, \end{aligned} \label{eq:BFSI-int-lin-discret} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{uu}= & \int_{\partial\Omega^{b}}\alpha_{f}N_{a}\frac{\varphi^{s}}{\varphi^{f}}\left(\boldsymbol{\tau}\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\right)\otimes\grad N_{b}\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\alpha_{f}N_{a}\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\cdot\boldsymbol{\mathcal{C}}_{d}\cdot\grad N_{b}\cdot\left(-\frac{\varphi^{s}}{\varphi^{f^{2}}}\mathbf{D}^{w}+\mathbf{M}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\alpha_{f}N_{a}\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(2\frac{\varphi^{s}}{\varphi^{f^{3}}}\grad\varphi^{f}\otimes\grad N_{b}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\alpha_{f}N_{a}\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\left(\frac{1}{\varphi^{f^{2}}}\grad N_{b}\otimes\grad\varphi^{f}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}\alpha_{f}N_{a}\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\cdot\boldsymbol{\mathcal{C}}^{\tau}\cdot\mathbf{w}\cdot\frac{\varphi^{s}}{\varphi^{f^{2}}}\left(-\frac{1}{J}\grad J\otimes\grad N_{b}+\grad\grad N_{b}\right)\thinspace d\eta^{1}d\eta^{2}\\ & +\int_{\partial\Omega^{b}}-\alpha_{f}N_{a}\left(-p\mathbf{I}+\boldsymbol{\tau}\right)\cdot\boldsymbol{\mathcal{A}}\left\{ -\mathbf{g}_{2}\frac{\partial N_{b}}{\partial\eta^{1}}+\mathbf{g}_{1}\frac{\partial N_{b}}{\partial\eta^{2}}\right\} \thinspace d\eta^{1}d\eta^{2}\,,\\ \mathbf{K}_{ab}^{uw}= & \int_{\partial\Omega^{b}}-\alpha_{f}N_{a}\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\cdot\boldsymbol{\mathcal{C}}_{d}\cdot\frac{1}{\varphi^{f}}\left(\grad N_{b}-\frac{1}{\varphi^{f}}N_{b}\grad\varphi^{f}\right)\thinspace d\eta^{1}d\eta^{2}\,,\\ \mathbf{k}_{ab}^{uJ}= & \int_{\partial\Omega^{b}}-\alpha_{f}N_{a}N_{b}\left(-p^{\prime}\mathbf{I}+\boldsymbol{\tau}_{J}^{\prime}\right)\cdot\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\thinspace d\eta^{1}d\eta^{2}\,. \end{aligned} \label{eq:BFSI-int-lin-discret-parts} \end{equation} \]
Note that \(\boldsymbol{\mathcal{A}}\left\{ \cdot\right\}\) represents the skew-symmetric tensor form of the bracketed vector.