3.2 Weak formulation for biphasic materials

A weak form of the statement conservation of linear momemtum for the quasi-static case is obtained by using Eqs.(2.7-2) and (2.7-4):

\[ \begin{equation} \delta W=\int_{b}\left[\delta\mathbf{v}^{s}\cdot\left(\divg\boldsymbol{\sigma}+\rho\mathbf{b}\right)+\delta p\,\divg\left(\mathbf{v}^{s}+\mathbf{w}\right)\right]dv=0\,,\label{eq193} \end{equation} \]

where \(b\) is the domain of interest defined on the solid matrix, \(\delta\mathbf{v}^{s}\) is a virtual velocity of the solid and \(\delta p\) is a virtual pressure of the fluid ¹. \(dv\) is an elemental volume of \(b\). Using the divergence theorem, this expression may be rearranged as

\[ \begin{equation} \begin{aligned}\delta W & =\int_{\partial b}\delta\mathbf{v}^{s}\cdot\mathbf{t}\,da+\int_{\partial b}\delta p\,w_{n}\,da+\int_{b}\delta\mathbf{v}^{s}\cdot\rho\mathbf{b}\,dv\\ & -\int_{b}\boldsymbol{\sigma}:\grad\delta\mathbf{v}^{s}\,dv-\int_{b}\left(\mathbf{w}\cdot\grad\delta p-\delta p\,\divg\mathbf{v}^{s}\right)\,dv \end{aligned} \,,\label{eq194} \end{equation} \]

where \(\delta\mathbf{d}^{s}=\left(\grad\delta\mathbf{v}^{s}+\grad^{T}\delta\mathbf{v}^{s}\right)/2\) is the virtual rate of deformation tensor, \(\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n}\) is the total traction on the surface \(\partial b\), and \(w_{n}=\mathbf{w}\cdot\mathbf{n}\) is the component of the fluid flux normal to \(\partial b\), with \(\mathbf{n}\) representing the unit outward normal to \(\partial b\). \(da\) represents an elemental area of \(\partial b\). In this type of problem, essential boundary conditions are prescribed for \(\mathbf{u}\) and \(p\), and natural boundary conditions are prescribed for \(\mathbf{t}\) and \(w_{n}\). In the expression of eq.\eqref{eq194}, \(\delta W\left(\boldsymbol{\chi}^{s},p,\delta\mathbf{v}^{s},\delta p\right)\) represents the virtual work.

Linearization

Since the system of equations in eq.\eqref{eq194} is highly nonlinear, its solution requires an iterative scheme such as Newton's method. This requires the linearization of \(\delta W\) at some trial solution \(\left(\boldsymbol{\chi}_{k}^{s},p_{k}\right)\), along an increment \(\Delta\mathbf{u}\) in \(\boldsymbol{\chi}^{s}\) and an increment \(\Delta p\) in \(p\),

\[ \begin{equation} \delta W+D\delta W\left[\Delta\mathbf{u}\right]+D\delta W\left[\Delta p\right]=0,\label{eq195} \end{equation} \]

where \(Df\left[\Delta q\right]\) represents the directional derivative of \(f\) along \(\Delta q\). For convenience, the virtual work may be separated into its internal and external parts,

\[ \begin{equation} \delta W=\delta W_{\text{ext}}-\delta W_{\text{int}},\label{eq196} \end{equation} \]

where

\[ \begin{equation} \delta W_{\mbox{int}}=\int_{b}\boldsymbol{\sigma}:\delta\mathbf{d}^{s}\,dv+\int_{b}\left(\mathbf{w}\cdot\grad\delta p-\delta p\,\frac{\dot{J}}{J}\right)\,dv,\label{eq197} \end{equation} \]

where we have substituted \(\divg\mathbf{v}^{s}=\dot{J}/J\), and

\[ \begin{equation} \delta W_{\text{ext}}=\int_{\partial b}\delta\mathbf{v}^{s}\cdot\mathbf{t}\,da+\int_{\partial b}\delta p\,w_{n}\,da+\int_{b}\delta\mathbf{v}^{s}\cdot\rho\mathbf{b}\,dv\,.\label{eq198} \end{equation} \]

The evaluation of the directional derivatives can be performed following a standard approach ². In particular, a backward difference scheme is used to evaluate \(\dot{J}\approx\left(J-J^{-\Delta t}\right)/\Delta t\), where \(J^{-\Delta t}\) is the value of \(J\) at the previous time step. For the internal part of the virtual work, the directional derivative along \(\Delta\mathbf{u}\) yields

\[ \begin{equation} \begin{aligned}D\delta W_{\text{int}}\left[\Delta\mathbf{u}\right] & =\int_{b}\delta\mathbf{d}^{s}:\boldsymbol{\mathcal{C}}:\Delta\boldsymbol{\varepsilon}\,dv+\int_{b}\boldsymbol{\sigma}:\left(\grad^{T}\Delta\mathbf{u}\cdot\grad\delta\mathbf{v}^{s}\right)\,dv\\ & -\int_{b}\frac{\delta p}{\Delta t}\divg\Delta\mathbf{u}\,dv\\ & -\int_{b}\grad\delta p\cdot\left(\boldsymbol{\mathcal{K}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\left(\grad p-\rho_{T}^{w}\mathbf{b}^{w}\right)\,dv\\ & +\int_{b}\grad\delta p\cdot\mathbf{k}\cdot\rho_{T}^{w}\left(\grad^{T}\Delta\mathbf{u}\cdot\mathbf{b}^{w}+\grad\mathbf{b}^{w}\cdot\Delta\mathbf{u}\right)dv\,, \end{aligned} \label{eq199} \end{equation} \]

where \(\boldsymbol{\mathcal{C}}\) is the fourth-order spatial elasticity tensor for the mixture and \(\Delta\boldsymbol{\varepsilon}=\left(\grad\Delta\mathbf{u}+\grad^{T}\Delta\mathbf{u}\right)/2\). Based on the relation of eq.(2.7-3), the spatial elasticity tensor may also be expanded as

\[ \begin{equation} \boldsymbol{\mathcal{C}}=\boldsymbol{\mathcal{C}}^{e}+p\left(-\mathbf{I}\otimes\mathbf{I}+2\mathbf{I}\,\overline{\underline{\otimes}}\,\mathbf{I}\right),\label{eq200} \end{equation} \]

where \(\boldsymbol{\mathcal{C}}^{e}\) is the spatial elasticity tensor for the solid matrix ³. It is related to the material elasticity tensor \(\boldsymbol{\mathbb{C}}^{e}\) via

\[ \begin{equation} \boldsymbol{\mathcal{C}}^{e}=J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):\boldsymbol{\mathbb{C}}^{e}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right),\label{eq201} \end{equation} \]

where \(\mathbf{F}\) is the deformation gradient of the solid matrix, \(\boldsymbol{\mathbb{C}}^{e}=\partial\mathbf{S}^{e}/\partial\mathbf{E}\) where \(\mathbf{E}\) is the Lagrangian strain tensor and \(\mathbf{S}^{e}\) is the second Piola-Kirchhoff stress tensor, related to the Cauchy stress tensor via \(\boldsymbol{\sigma}^{e}=J^{-1}\mathbf{F}\cdot\mathbf{S}^{e}\cdot\mathbf{F}^{T}\).

Similarly, \(\boldsymbol{\mathcal{K}}\) is a fourth-order tensor that represents the spatial measure of the rate of change of permeability with strain. It is related to its material frame equivalent \(\boldsymbol{\mathbb{K}}\) via

\[ \begin{equation} \boldsymbol{\mathcal{K}}=J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):\boldsymbol{\mathbb{K}}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right),\label{eq202} \end{equation} \]

where \(\boldsymbol{\mathbb{K}}=\partial\mathbf{K}/\partial\mathbf{E}\) and \(\mathbf{K}\) is the permeability tensor in the material frame, such that \(\mathbf{k}=J^{-1}\mathbf{F}\cdot\mathbf{K}\cdot\mathbf{F}^{T}\). Since \(\mathbf{K}\) and \(\mathbf{E}\) are symmetric tensors, it follows that \(\boldsymbol{\mathcal{K}}\) and \(\boldsymbol{\mathbb{K}}\) exhibit two minor symmetries (e.g., \(\mathcal{K}_{jikl}=\mathcal{K}_{ijkl}\) and \(\mathcal{K}_{ijlk}=\mathcal{K}_{ijkl})\). However, unlike the elasticity tensor, it is not necessary that these tensors exhibit major symmetry (e.g., \(\mathcal{K}_{klij}\ne\mathcal{K}_{ijkl}\) in general).

The directional derivative of \(\delta W_{\mbox{int}}\) along \(\Delta p\) is given by

\[ \begin{equation} D\delta W_{\text{int}}\left[\Delta p\right]=-\int_{b}\grad\delta p\cdot\mathbf{k}\cdot\grad\Delta p\,dv-\int_{b}\Delta p\,\divg\delta\mathbf{v}^{s}\,dv\,.\label{eq203} \end{equation} \]

Note that letting \(p=0\) and \(\delta p=0\) in the above equations recovers the virtual work relations for nonlinear elasticity of compressible solids. The resulting simplified equation emerging from eq.\eqref{eq199} is symmetric to interchanges of \(\Delta\mathbf{u}\) and \(\delta\mathbf{v}^{s}\), producing a symmetric stiffness matrix in the finite element formulation ². However, the general relations of Eqs.\eqref{eq199} and \eqref{eq203} do not exhibit symmetry to interchanges of \(\left(\Delta\mathbf{u},\Delta p\right)\) and \(\left(\delta\mathbf{v}^{s},\delta p\right)\), implying that the finite element stiffness matrix for a solid-fluid mixture is not symmetric under general conditions.

The directional derivatives of the external virtual work \(\delta W_{\text{ext}}\) depend on the type of boundary conditions being considered. For a prescribed total normal traction \(t_{n}\), where \(\mathbf{t}=t_{n}\mathbf{n}\),

\[ \begin{equation} \delta W_{\text{ext}}^{t}=\int_{\partial b}\delta\mathbf{v}^{s}\cdot t_{n}\mathbf{n}\,da,\label{eq204} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}D\delta W_{\text{ext}}^{t}\left[\Delta\mathbf{u}\right] & =\int_{\partial b}\delta\mathbf{v}^{s}\cdot t_{n}\left(\mathbf{g}_{1}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{2}}-\mathbf{g}_{2}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{1}}\right)\frac{da}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,,\\ D\delta W_{\text{ext}}^{t}\left[\Delta p\right] & =0, \end{aligned} \label{eq205} \end{equation} \]

where

\[ \begin{equation} \mathbf{g}_{\alpha}=\frac{\partial\mathbf{x}}{\partial\eta^{\alpha}},\quad\alpha=1,2\label{eq206} \end{equation} \]

are covariant basis (tangent) vectors on \(\partial b\), such that

\[ \begin{equation} \mathbf{n}=\frac{\mathbf{g}_{1}\times\mathbf{g}_{2}}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,.\label{eq207} \end{equation} \]

For a prescribed normal effective traction \(t_{n}^{e}\), where \(\mathbf{t}=\left(-p+t_{n}^{e}\right)\mathbf{n}\) and \(p\) is not prescribed, then

\[ \begin{equation} \delta W_{\text{ext}}^{e}=\int_{\partial b}\delta\mathbf{v}^{s}\cdot\left(-p+t_{n}^{e}\right)\mathbf{n}\,da\,,\label{eq208} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}D\delta W_{\text{ext}}^{e}\left[\Delta\mathbf{u}\right] & =\int_{\partial b}\delta\mathbf{v}^{s}\cdot\left(-p+t_{n}^{e}\right)\left(\mathbf{g}_{1}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{2}}-\mathbf{g}_{2}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{1}}\right)\frac{da}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,,\\ D\delta W_{\text{ext}}^{e}\left[\Delta p\right] & =-\int_{\partial b}\delta\mathbf{v}^{s}\cdot\Delta p\,\mathbf{n}\,da\,. \end{aligned} \label{eq209} \end{equation} \]

For a prescribed normal fluid flux \(w_{n}=\mathbf{w}\cdot\mathbf{n}\),

\[ \begin{equation} \delta W_{\text{ext}}^{w}=\int_{\partial b}\delta p\,w_{n}\,da\,,\label{eq210} \end{equation} \]

and

\[ \begin{equation} \begin{aligned}D\delta W_{\text{ext}}^{w}\left[\Delta\mathbf{u}\right] & =\int_{\partial b}\delta p\,w_{n}\,\mathbf{n}\cdot\left(\mathbf{g}_{1}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{2}}-\mathbf{g}_{2}\times\frac{\partial\Delta\mathbf{u}}{\partial\eta^{1}}\right)\frac{da}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,,\\ D\delta W_{\text{ext}}^{w}\left[\Delta p\right] & =0\,. \end{aligned} \label{eq211} \end{equation} \]

Finally, for a prescribed external body force, recognizing that \(\rho\mathbf{b}=\rho^{s}\mathbf{b}^{s}+\rho^{w}\mathbf{b}^{w}\) and assuming that the body forces \(\mathbf{b}^{s}\) and \(\mathbf{b}^{w}\) do not depend on \(p\),

\[ \begin{equation} \begin{aligned}D\left(\delta W_{\text{ext}}^{b}\right)\left[\Delta\mathbf{u}\right] & =\int_{b}\delta\mathbf{v}^{s}\cdot\left[\left(\rho^{s}grad\mathbf{b}^{s}+\rho^{w}grad\mathbf{b}^{w}\right)\cdot\Delta\mathbf{u}+\left(\divg\Delta\mathbf{u}\right)\rho_{T}^{w}\mathbf{b}^{w}\right]\,dv\,,\\ D\left(\delta W_{\text{ext}}^{b}\right)\left[\Delta p\right] & =0\,. \end{aligned} \label{eq212} \end{equation} \]

Stabilization

Finite element models of deformable porous media are known to exhibit oscillations in fluid pressure caused by relatively coarse meshes near free-draining boundaries, where a boundary layer in fluid pressure normally forms during the early time response to loading. These spurious oscillations typically occur when the mesh is not able to resolve this boundary layer. Stabilization methods serve as an alternative to mesh refinement, when the latter is not feasible or practical. In FEBio we have implemented a stabilization method based on the work of Aguilar et al. ⁴. Basically, this method proposes that the fluid flux \(\mathbf{w}\) be evaluated from

\[ \begin{equation} \mathbf{w}=-\mathbf{k}\cdot\grad\left(p+\tau\dot{p}\right)\label{eq:stab-fluid-flux} \end{equation} \]

instead of eq.(2.7-6), where \(\tau\) represents the stabilization parameter. A representative value for \(\tau\) may be evaluated from

\[ \begin{equation} \tau\approx\frac{h^{2}}{E\cdot k}\,,\label{eq:stab-tau} \end{equation} \]

where \(h\) is the element thickness on the free-draining boundary, \(E\) is a representative measure of the solid matrix elastic modulus and \(k\) is a representative measure of the hydraulic permeability tensor \(\mathbf{k}\). The contribution of this form of \(\mathbf{w}\) to the virtual work expression in eq.\eqref{eq194} is

\[ \begin{equation} \begin{aligned}\delta G & =\int_{b}\mathbf{w}\cdot\grad\delta p\,dv\\ & =-\int_{b}\grad\delta p\cdot\mathbf{k}\cdot\grad\left(p+\tau\dot{p}\right)\,dv \end{aligned} \,.\label{eq:stab-virtual-work} \end{equation} \]

The linearizations of \(\delta G\) are then given by

\[ \begin{equation} \begin{aligned}D\delta G\left[\Delta\mathbf{u}\right] & =-\int_{b}\grad\delta p\cdot\left(\boldsymbol{\mathcal{K}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\grad\left(p+\tau\dot{p}\right)\,dv\\ D\delta G\left[\Delta p\right] & =-\int_{b}\left(1+\frac{\tau}{\Delta t}\right)\grad\delta p\cdot\mathbf{k}\cdot\grad\Delta p\,dv \end{aligned} \,.\label{eq:stab-linearization} \end{equation} \]

Discretization

Let

\[ \begin{equation} \begin{aligned}\delta\mathbf{v}^{s} & =\sum\limits_{a=1}^{m}N_{a}\delta\mathbf{v}_{a}\,, & \delta p & =\sum\limits_{a=1}^{m}N_{a}\delta p_{a}\,,\\ \Delta\mathbf{u} & =\sum\limits_{b=1}^{m}N_{b}\Delta\mathbf{u}_{b}\,, & \Delta p & =\sum\limits_{b=1}^{m}N_{b}\Delta p_{b}\,, \end{aligned} \label{eq213} \end{equation} \]

where \(N_{a}\) represents the interpolation functions over an element, \(\delta\mathbf{v}_{a},\delta p_{a},\Delta\mathbf{u}_{b}\mbox{\thinspace and\thinspace}\Delta p_{b}\) respectively represent nodal values of \(\delta\mathbf{v}^{s}\), \(\delta p\), \(\Delta\mathbf{u}\) and \(\Delta p\), and \(m\) is the number of nodes in an element. Then the discretized form of \(\delta W_{int}\) in eq.\eqref{eq197} may be written as

\[ \begin{equation} \delta W_{\text{int}}=\sum\limits_{e=1}^{n_{e}}\sum\limits_{k=1}^{n_{\text{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\left[\begin{array}{cc} \delta\mathbf{v}_{a} & \delta p_{a}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{r}_{a}^{u}\\ r_{a}^{p} \end{array}\right]\,,\label{eq214} \end{equation} \]

where \(n_{e}\) is the number of elements in \(b\), \(n_{\mbox{int}}^{\left(e\right)}\) is the number of integration points in the \(e-\)th element, \(W_{k}\) is the quadrature weight associated with the \(k-\)th integration point, and \(J_{\eta}\) is the Jacobian of the transformation from the spatial frame to the parametric space of the element. In the above expression,

\[ \begin{equation} \mathbf{r}_{a}^{u}=\boldsymbol{\sigma}\cdot\nabla N_{a}\,,\quad r_{a}^{p}=\mathbf{w}\cdot\nabla N_{a}-N_{a}\mbox{div}\mathbf{v}^{s}\,,\label{eq215} \end{equation} \]

and it is understood that \(J_{\eta}\), \(\mathbf{r}_{a}^{u}\) and \(r_{a}^{p}\) are evaluated at the parametric coordinates of the \(k-\)th integration point.

Similarly, the discretized form of \(D\delta W_{\mbox{int}}\) in Eqs.\eqref{eq199} and \eqref{eq203} may be written as

\[ \begin{equation} -D\delta W_{\text{int}}=\sum\limits_{e=1}^{n_{e}}\sum\limits_{k=1}^{n_{\text{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\left[\begin{array}{cc} \delta\mathbf{v}_{a} & \delta p_{a}\end{array}\right]\cdot\sum\limits_{b=1}^{m}\left[\begin{array}{cc} \mathbf{K}_{ab}^{uu} & \mathbf{k}_{ab}^{up}\\ \mathbf{k}_{ab}^{pu} & k_{ab}^{pp} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta p_{b} \end{array}\right]\,,\label{eq216} \end{equation} \]

where

\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{uu} & =\nabla N_{a}\cdot\boldsymbol{\mathcal{C}}\cdot\nabla N_{b}+\left(\nabla N_{a}\cdot\boldsymbol{\sigma}\cdot\nabla N_{b}\right)\mathbf{I}\\ & -N_{a}\left[N_{b}\left(\rho^{s}\nabla\mathbf{b}^{s}+\rho^{w}\nabla\mathbf{b}^{w}\right)+\rho_{T}^{w}\mathbf{b}^{w}\otimes\nabla N_{b}\right]\,,\\ \mathbf{k}_{ab}^{up} & =-N_{b}\nabla N_{a}\,,\\ \mathbf{k}_{ab}^{pu} & =-\left(\nabla N_{a}\cdot\boldsymbol{\mathcal{K}}\cdot\nabla N_{b}\right)\cdot\left(\nabla p-\rho_{T}^{w}\mathbf{b}^{w}\right)-\frac{1}{\Delta t}N_{a}\cdot\nabla N_{b}\\ & +\rho_{T}^{w}\left(\mathbf{b}^{w}\otimes\nabla N_{b}+N_{b}\nabla^{T}\mathbf{b}^{w}\right)\cdot\mathbf{k}\cdot\nabla N_{a}\,,\\ k_{ab}^{pp} & =-\nabla N_{a}\cdot\mathbf{k}\cdot\nabla N_{b}, \end{aligned} \label{eq217} \end{equation} \]

and \(\Delta t\) is a discrete increment in time. In a numerical implementation, it has been found that evaluating \(\divg\mathbf{v}^{s}\) from \(\dot{J}/J\), where \(J=\det\mathbf{F}\), yields more accurate solutions than evaluating it from the trace of \(\grad\mathbf{v}^{s}\) ⁵. Contributions from the stabilization method presented in Section \eqref{subsec:Biphasic-Stabilization} may be similarly evaluated.

For the various types of contributions to the external virtual work, a similar discretization produces

\[ \begin{equation} \delta W_{\text{ext}}=\sum\limits_{e=1}^{n_{e}}\sum\limits_{k=1}^{n_{\text{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\left[\begin{array}{cc} \delta\mathbf{v}_{a} & \delta p_{a}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{r}_{a}^{u}\\ r_{a}^{p} \end{array}\right]\,,\label{eq218} \end{equation} \]

and

\[ \begin{equation} -D\delta W_{\text{ext}}=\sum\limits_{e=1}^{n_{e}}\sum\limits_{k=1}^{n_{\text{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\left[\begin{array}{cc} \delta\mathbf{v}_{a} & \delta p_{a}\end{array}\right]\cdot\sum\limits_{b=1}^{m}\left[\begin{array}{cc} \mathbf{K}_{ab}^{uu} & \mathbf{k}_{ab}^{up}\\ \mathbf{k}_{ab}^{pu} & k_{ab}^{pp} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta p_{b} \end{array}\right]\,,\label{eq219} \end{equation} \]

where

\[ \begin{equation} J_{\eta}=\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\,.\label{eq220} \end{equation} \]

In this case, \(m\) represents the number of nodes on an element face. For a prescribed normal traction \(t_{n}\) as given in \eqref{eq204}-\eqref{eq205},

\[ \begin{equation} \begin{aligned}\mathbf{r}_{a}^{u} & =t_{n}N_{a}\mathbf{n}\,, & r_{a}^{u} & =0\,,\\ \mathbf{K}_{ab}^{uu} & =t_{n}N_{a}\frac{1}{J_{\eta}}\boldsymbol{\mathbf{\mathcal{A}}}\left\{ \frac{\partial N_{b}}{\partial\eta^{1}}\mathbf{g}_{2}-\frac{\partial N_{b}}{\partial\eta^{2}}\mathbf{g}_{1}\right\} \,, & \mathbf{k}_{ab}^{up} & =\mathbf{0}\,,\\ \mathbf{k}_{ab}^{pu} & =\mathbf{0}\,, & k_{ab}^{pp} & =0\,, \end{aligned} \label{eq221} \end{equation} \]

where \(\boldsymbol{\mathbf{\mathcal{A}}}\left\{ \mathbf{v}\right\} =-\boldsymbol{\mathbf{\mathcal{E}}}\cdot\mathbf{v}\) is the skew-symmetric tensor whose dual vector is \(\mathbf{v}\) and \(\boldsymbol{\mathbf{\mathcal{E}}}\) is the third-order permutation pseudo-tensor. For a prescribed traction \(t_{n}^{e}\) as given in \eqref{eq208}-\eqref{eq209},

\[ \begin{equation} \begin{aligned}\mathbf{r}_{a}^{u} & =\left(-p+t_{n}^{e}\right)N_{a}\mathbf{n}\,, & r_{a}^{u} & =0\,,\\ \mathbf{K}_{ab}^{uu} & =\left(-p+t_{n}^{e}\right)N_{a}\frac{1}{J_{\eta}}\boldsymbol{\mathbf{\mathcal{A}}}\left\{ \frac{\partial N_{b}}{\partial\eta^{1}}\mathbf{g}_{2}-\frac{\partial N_{b}}{\partial\eta^{2}}\mathbf{g}_{1}\right\} \,, & \mathbf{k}_{ab}^{up} & =\mathbf{0}\,,\\ \mathbf{k}_{ab}^{pu} & =\mathbf{0}\,, & k_{ab}^{pp} & =0\,. \end{aligned} \label{eq222} \end{equation} \]

For a prescribed normal fluid flux \(w_{n}\) as given in \eqref{eq210}-\eqref{eq211},

\[ \begin{equation} \begin{aligned}\mathbf{r}_{a}^{u} & =\mathbf{0}\,, & r_{a}^{u} & =w_{n}N_{a}\,,\\ \mathbf{K}_{ab}^{uu} & =\mathbf{0}\,, & \mathbf{k}_{ab}^{up} & =\mathbf{0}\,,\\ \mathbf{k}_{ab}^{pu} & =w_{n}N_{a}\frac{1}{J_{\eta}}\mathbf{n}\times\left(\frac{\partial N_{b}}{\partial\eta^{1}}\mathbf{g}_{2}-\frac{\partial N_{b}}{\partial\eta^{2}}\mathbf{g}_{1}\right)\,, & k_{ab}^{pp} & =0\,. \end{aligned} \label{eq223} \end{equation} \]

---

1. Un, K.; Spilker, R. L.. "A penetration-based finite element method for hyperelastic 3D biphasic tissues in contact. Part II: finite element simulations." J Biomech Eng, vol. 128, pp. 934-42 (2006). ↩

2. Bonet, Javier; Wood, Richard D.. "Nonlinear continuum mechanics for finite element analysis." Cambridge University Press (1997). ↩↩

3. Curnier, A.; Qi-Chang, He; Zysset, P.. "Conewise linear elastic materials." J Elasticity, vol. 37, pp. 1-38 (1994). ↩

4. Aguilar, G; Gaspar, F; Lisbona, F; Rodrigo, C24558811158. "Numerical stabilization of Biot's consolidation model by a perturbation on the flow equation." International journal for numerical methods in engineering, vol. 75, pp. 1282--1300 (2008). ↩

5. Ateshian, Gerard A; Ellis, Benjamin J; Weiss, Jeffrey A. "Equivalence between short-time biphasic and incompressible elastic material responses." J Biomech Eng, vol. 129, pp. 405-12 (2007). ↩