The virtual work integral for this problem is given by
\[ \begin{equation} \begin{aligned}\delta W & =\int_{b}\delta\mathbf{v}\cdot\divg\boldsymbol{\sigma}\,dv+\int_{b}\delta\tilde{p}\,\divg\left(\mathbf{v}^{s}+\mathbf{w}\right)\,dv\\ & +\int_{b}\delta\tilde{c}\left[\frac{\partial\left(\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}+\divg\left(\mathbf{j}+\phi^{w}\tilde{\kappa}\tilde{c}\,\mathbf{v}^{s}\right)\right]\,dv\,, \end{aligned} \label{eq224} \end{equation} \]
where \(\delta\mathbf{v}\) is the virtual velocity of the solid, \(\delta\tilde{p}\) is the virtual effective fluid pressure, and \(\delta\tilde{c}\) is the virtual molar energy of the solute. \(b\) represents the mixture domain in the spatial frame and \(dv\) is an elemental mixture volume in \(b\). In the last integral of \(\delta W\), note that
\[ \begin{equation} \frac{\partial\left(\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}+\divg\left(\varphi^{w}\tilde{\kappa}\tilde{c}\,\mathbf{v}^{s}\right)=\frac{1}{J}\frac{D^{s}}{Dt}\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)\,,\label{eq225} \end{equation} \]
where \(D^{s}f/Dt\equiv\partial f/\partial t+\mathbf{v}^{s}\cdot\grad f\) is the material time derivative of a scalar function \(f\) in the spatial frame, following the solid. Similarly, note that \(\divg\mathbf{v}^{s}=J^{-1}\left(D^{s}J/Dt\right)\). Using the divergence theorem, the virtual work integral may be separated into internal and external contributions, \(\delta W=\delta W_{\text{ext}}-\delta W_{\text{int}}\), where
\[ \begin{equation} \begin{aligned}\delta W_{\text{int}} & =\int_{b}\boldsymbol{\sigma}:\delta\mathbf{d}^{s}\,dv+\int_{b}\left(\mathbf{w}\cdot\grad\delta\tilde{p}-\frac{\delta\tilde{p}}{J}\frac{D^{s}J}{Dt}\right)\,dv\\ & +\int_{b}\left[\mathbf{j}\cdot\grad\delta\tilde{c}-\frac{\delta\tilde{c}}{J}\frac{D^{s}}{Dt}\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)\right]dv\,,\\ \delta W_{\text{ext}} & =\int_{\partial b}\left(\delta\mathbf{v}\cdot\mathbf{t}+\delta\tilde{p}\,w_{n}+\delta\tilde{c}\,j_{n}\right)\,da\,, \end{aligned} \label{eq226} \end{equation} \]
with \(\delta W_{\text{ext}}\) being evaluated on the domain's boundary surface \(\partial b\). In the first expression \(\delta\mathbf{d}^{s}=\left(\grad\delta\mathbf{v}+\grad^{T}\delta\mathbf{v}\right)/2\) represents the virtual solid rate of deformation.
To solve this nonlinear system using an iterative Newton scheme, the virtual work must be linearized at trial solutions, along increments in \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}\),
\[ \begin{equation} \delta W+D\delta W\left[\Delta\mathbf{u}\right]+D\delta W\left[\Delta\tilde{p}\right]+D\delta W\left[\Delta\tilde{c}\right]\approx0\,,\label{eq227} \end{equation} \]
where, for any function \(f\left(q\right)\), \(Df\left[\Delta q\right]\) represents the directional derivative of \(f\) along \(\Delta q\) . To operate the directional derivative on the integrand of \(\delta W_{\text{int}}\), it is first necessary to convert the integrals from the spatial to the material domain :
\[ \begin{equation} \delta W_{\text{int}}=\int_{B}\mathbf{S}:\delta\mathbf{\dot{E}}\,dV+\int_{B}\left(\mathbf{W}\cdot\grad\delta\tilde{p}-\delta\tilde{p}\frac{\partial J}{\partial t}\right)\,dV+\int_{B}\left[\mathbf{J}\cdot\grad\delta\tilde{c}-\delta\tilde{c}\frac{\partial}{\partial t}\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)\right]dV\,,\label{eq228} \end{equation} \]
where \(B\) represents the mixture domain in the material frame, \(dV\) is an elemental mixture volume in \(B\), and
\[ \begin{equation} \begin{aligned}\mathbf{S} & =J\,\mathbf{F}^{-1}\cdot\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\,,\\ \delta\mathbf{\dot{E}} & =\mathbf{F}^{T}\cdot\delta\mathbf{d}^{s}\cdot\mathbf{F}\,,\\ \mathbf{W} & =J\,\mathbf{F}^{-1}\cdot\mathbf{w}\,,\\ \mathbf{J} & =J\,\mathbf{F}^{-1}\cdot\mathbf{j}\,. \end{aligned} \label{eq229} \end{equation} \]
The second Piola-Kirchhoff stress tensor \(\mathbf{S}\), and material flux vectors \(\mathbf{W}\) and \(\mathbf{J}\), are respectively related to \(\boldsymbol{\sigma}\), \(\mathbf{w}\) and \(\mathbf{j}\) by the Piola transformations for tensors and vectors . Substituting \eqref{eq229} into (2.8-13) produces
\[ \begin{equation} \begin{aligned}\mathbf{W} & =-\mathbf{\tilde{K}}\cdot\left(\grad\tilde{p}+R\theta\frac{\tilde{\kappa}}{d_{0}}J^{-1}\mathbf{C}\cdot\mathbf{D}\cdot\grad\tilde{c}\right)\,,\\ \mathbf{J} & =\tilde{\kappa}\mathbf{D}\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}J^{-1}\mathbf{C}\cdot\mathbf{W}\right)\,, \end{aligned} \label{eq230} \end{equation} \]
where \(\mathbf{\tilde{K}}\) and \(\mathbf{D}\) are the material representations of the permeability and diffusivity tensors, related to \(\mathbf{\tilde{k}}\) and \(\mathbf{d}\) via the Piola transformation,
\[ \begin{equation} \begin{aligned}\mathbf{\tilde{K}} & =J\,\mathbf{F}^{-1}\cdot\mathbf{\tilde{k}}\cdot\mathbf{F}^{-T}\,,\\ \mathbf{D} & =J\,\mathbf{F}^{-1}\cdot\mathbf{d}\cdot\mathbf{F}^{-T}\,. \end{aligned} \label{eq231} \end{equation} \]
The linearization of \(\delta W_{\text{int}}\) is rather involved and a summary of the resulting lengthy expressions is provided below. In consideration of the dearth of experimental data relating \(\tilde{\kappa}\) and \(\Phi\) to the complete state of solid matrix strain (such as \(\mathbf{C})\), this implementation assumes that the dependence of these functions on the strain is restricted to a dependence on the volume ratio \(J=\left(\det\mathbf{C}\right)^{1/2}\). Furthermore, it is assumed that the free solution diffusivity \(d_{0}\) is independent of the strain.
The linearization of \(\delta W_{\text{ext}}\) is described in Section Linearization of External Virtual Work. Following the linearization procedure, the resulting expressions may be discretized by nodally interpolating \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}\) over finite elements, producing a set of equations in matrix form, as described in Section Linearization of External Virtual Work.
The formulation presented in this study is implemented in FEBio by introducing an additional module dedicated to solute transport in deformable porous media. Classes are implemented to describe material functions for \(\boldsymbol{\sigma}^{e}\), \(\mathbf{k}\), \(\mathbf{d}\) (and \(d_{0})\), \(\tilde{\kappa}\) and \(\Phi\), which allow the formulation of any desired constitutive relation for these functions of \(\mathbf{C}\) and \(\tilde{c}\), along with corresponding derivatives of these functions with respect to \(\mathbf{C}\) and \(\tilde{c}\). The implementation accepts essential boundary conditions on \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}\), or natural boundary conditions on \(\mathbf{t}\), \(w_{n}\) and \(j_{n}\); initial conditions may also be specified for \(\tilde{p}\) and \(\tilde{c}\). Analysis results for pressure and concentration may be displayed either as \(\tilde{p}\) and \(\tilde{c}\), or as \(p\) and \(c\) by inverting the relations of (2.8-11).
Linearization of Internal Virtual Work
The virtual work integral \(\delta W_{\text{int}}\) in \eqref{eq228} may be linearized term by term along increments in \(\Delta\mathbf{u}\), \(\Delta\tilde{p}\) and \(\Delta\tilde{c}\) using the general form
\[ \begin{equation} D\left(\int_{B}F\,dV\right)\left[\Delta q\right]=\int_{B}DF\left[\Delta q\right]\,dV=\int_{b}f\,dv\,.\label{eq232} \end{equation} \]
For notational simplicity, the integral sign is omitted and the linearization of each term is presented in the form \(DF\left[\Delta q\right]\,dV=f\,dv\).
Linearization along \(\Delta\mathbf{u}\)
The linearization of the first term in \(\delta W_{\text{int}}\) along \(\Delta\mathbf{u}\) yields
\[ \begin{equation} \left(\mathbf{S}:\delta\dot{\mathbf{E}}\right)\left[\Delta\mathbf{u}\right]\,dV=\left[\delta\mathbf{d}^{s}:\boldsymbol{\mathcal{C}}:\Delta\boldsymbol{\varepsilon}+\boldsymbol{\sigma}:\left(\grad^{T}\Delta\mathbf{u}\cdot\grad\delta\mathbf{v}\right)\right]\,dv\,,\label{eq233} \end{equation} \]
where \(\boldsymbol{\mathcal{C}}\) is the spatial elasticity tensor of the mixture,
\[ \begin{equation} \boldsymbol{\mathcal{C}}=\boldsymbol{\mathcal{C}}^{e}-\left(\tilde{p}+R\theta\,\Phi\tilde{\kappa}\tilde{c}\right)\left(\mathbf{I}\otimes\mathbf{I}-2\mathbf{I}\,\overline{\underline{\otimes}}\,\mathbf{I}\right)-R\theta\tilde{c}\,J\frac{\partial\left(\Phi\tilde{\kappa}\right)}{\partial J}\mathbf{I}\otimes\mathbf{I}\,,\label{eq234} \end{equation} \]
and \(\boldsymbol{\mathcal{C}}^{e}\) is the spatial elasticity tensor of the solid matrix,
\[ \begin{equation} \boldsymbol{\mathcal{C}}^{e}=J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):2\frac{\partial\mathbf{S}^{e}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right)\,.\label{eq235} \end{equation} \]
The linearization of the second term is
\[ \begin{equation} D\left(\mathbf{W}\cdot\Grad\delta\tilde{p}\right)\left[\Delta\mathbf{u}\right]\,dV=\grad\delta\tilde{p}\cdot\mathbf{w}_{u}^{\prime}\,dv\,,\label{eq236} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{u}^{\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{W}\left[\Delta\mathbf{u}\right]=-\left(\tilde{\boldsymbol{\mathcal{K}}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\left(\grad\tilde{p}+R\theta\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad\tilde{c}\right)\\ & -\frac{R\theta}{d_{0}}\tilde{\mathbf{k}}\cdot\left(J^{2}\frac{\partial\left(J^{-1}\tilde{\kappa}\right)}{\partial J}\left(\divg\Delta\mathbf{u}\right)\mathbf{I}+2\tilde{\kappa}\,\Delta\boldsymbol{\varepsilon}\right)\cdot\mathbf{d}\cdot\grad\tilde{c}-\tilde{\kappa}\frac{R\theta}{d_{0}}\tilde{\mathbf{k}}\cdot\left(\boldsymbol{\mathcal{D}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\grad\tilde{c}\,, \end{aligned} \label{eq237} \end{equation} \]
with
\[ \begin{equation} \begin{aligned}\tilde{\boldsymbol{\mathcal{K}}} & =J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):2\frac{\partial\tilde{\mathbf{K}}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right)\,,\\ \boldsymbol{\mathcal{D}} & =J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):2\frac{\partial\mathbf{D}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right)\,, \end{aligned} \label{eq238} \end{equation} \]
representing the spatial tangents, with respect to the strain, of the effective permeability and solute diffusivity, respectively. These fourth-order tensors exhibit minor symmetries but not major symmetry, as described recently . Since \(\tilde{\mathbf{K}}\) is given by substituting (2.8-13)\(_{3}\) into \eqref{eq231}\(_{1}\), the evaluation of \(\tilde{\boldsymbol{\mathcal{K}}}\) is rather involved and it can be shown that
\[ \begin{equation} \tilde{\boldsymbol{\mathcal{K}}}=2\left(\tilde{\mathbf{k}}\otimes\mathbf{I}-2\tilde{\mathbf{k}}\,\underline{\otimes}\,\mathbf{I}\right)-\left(\tilde{\mathbf{k}}\,\overline{\underline{\otimes}}\,\tilde{\mathbf{k}}\right):\boldsymbol{\mathcal{G}}\,,\label{eq239} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\boldsymbol{\mathcal{G}} & =2\left(\mathbf{k}^{-1}\otimes\mathbf{I}-2\mathbf{k}^{-1}\,\overline{\underline{\otimes}}\,\mathbf{I}\right)-\left(\mathbf{k}^{-1}\,\underline{\otimes}\,\mathbf{k}^{-1}\right):\boldsymbol{\mathcal{K}}\\ & +\frac{R\theta\tilde{c}}{d_{0}}J\frac{\partial}{\partial J}\left(\frac{\tilde{\kappa}}{\varphi^{w}}\right)\left(\mathbf{I}-\frac{\mathbf{d}}{d_{0}}\right)\otimes\mathbf{I}\\ & +\frac{R\theta\tilde{c}}{d_{0}}\frac{\tilde{\kappa}}{\varphi^{w}}\left(\mathbf{I}\otimes\mathbf{I}-2\mathbf{I}\,\underline{\otimes}\,\mathbf{I}-\frac{1}{d_{0}}\boldsymbol{\mathcal{D}}\right) \end{aligned} \label{eq240} \end{equation} \]
and
\[ \begin{equation} \boldsymbol{\mathcal{K}}=J^{-1}\left(\mathbf{F}\,\underline{\otimes}\,\mathbf{F}\right):2\frac{\partial\mathbf{K}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\,\underline{\otimes}\,\mathbf{F}^{T}\right)\,.\label{eq241} \end{equation} \]
The next term in \(\delta W_{\text{int}}\) linearizes to
\[ \begin{equation} -D\left(\delta\tilde{p}\frac{\partial J}{\partial t}\right)\left[\Delta\mathbf{u}\right]\,dV=-\delta\tilde{p}\frac{1}{\Delta t}\divg\Delta\mathbf{u}\,dv\,,\label{eq242} \end{equation} \]
where we used a backward difference scheme to approximate the time derivative,
\[ \begin{equation} \frac{\partial J}{\partial t}\approx\frac{1}{\Delta t}\left(J-J^{-\Delta t}\right)\,,\label{eq243} \end{equation} \]
and \(\Delta t\) represents the time increment relative to the previous time point. The next term is given by
\[ \begin{equation} D\left(\mathbf{J}\cdot\Grad\delta\tilde{c}\right)\left[\Delta\mathbf{u}\right]\,dV=\grad\delta\tilde{c}\cdot\mathbf{j}_{u}^{\prime}\,dv,\label{eq244} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{j}_{u}^{\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{J}\left[\Delta\mathbf{u}\right]=\left(J\frac{\partial\tilde{\kappa}}{\partial J}\left(\divg\Delta\mathbf{u}\right)\mathbf{d}+\tilde{\kappa}\boldsymbol{\mathcal{D}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}\right)\\ & +\tilde{\kappa}\mathbf{d}\cdot\left(-\varphi^{s}\left(\divg\Delta\mathbf{u}\right)\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\left(2\Delta\boldsymbol{\varepsilon}-\left(\divg\Delta\mathbf{u}\right)\mathbf{I}\right)\cdot\mathbf{w}\right)+\tilde{\kappa}\frac{\tilde{c}}{d_{0}}\mathbf{d}\cdot\mathbf{w}_{u}^{\prime} \end{aligned} \,.\label{eq245} \end{equation} \]
Using a backward difference scheme for the time derivative, the last term is
\[ \begin{equation} D\left(\delta\tilde{c}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}\right)\left[\Delta\mathbf{u}\right]\,dV=\frac{\delta\tilde{c}}{\Delta t}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}\right)}{\partial J}\tilde{c}\divg\Delta\mathbf{u}\,dv\,.\label{eq246} \end{equation} \]
Linearization along \(\Delta\tilde{p}\)
The linearization of the various terms in \(\delta W_{\text{int}}\) along \(\Delta\tilde{p}\) yields
\[ \begin{equation} D\left(\mathbf{S}:\delta\mathbf{\dot{E}}\right)\left[\Delta\tilde{p}\right]\,dV=-\Delta\tilde{p}\,\divg\delta\mathbf{v}\,dv\,,\label{eq247} \end{equation} \]
\[ \begin{equation} D\left(\mathbf{W}\cdot\Grad\delta\tilde{p}-\delta\tilde{p}\frac{\partial J}{\partial t}\right)\left[\Delta\tilde{p}\right]\,dV=-\grad\delta\tilde{p}\cdot\tilde{\mathbf{k}}\cdot\grad\Delta\tilde{p}\,dv\,,\label{eq248} \end{equation} \]
\[ \begin{equation} D\left(\mathbf{J}\cdot\Grad\delta\tilde{c}-\delta\tilde{c}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}\right)\left[\Delta\tilde{p}\right]\,dV=-\frac{\tilde{\kappa}\tilde{c}}{d_{0}}\grad\delta\tilde{c}\cdot\mathbf{d}\cdot\tilde{\mathbf{k}}\cdot\grad\Delta\tilde{p}\,dv\,.\label{eq249} \end{equation} \]
Linearization along \(\Delta\tilde{c}\)
The linearization of the first term in \(\delta W_{\text{int}}\) along \(\Delta\tilde{c}\) yields
\[ \begin{equation} D\left(\mathbf{S}:\delta\mathbf{\dot{E}}\right)\left[\Delta\tilde{c}\right]\,dV=\Delta\tilde{c}\left(\boldsymbol{\sigma}_{c}^{\prime}:\delta\mathbf{d}-R\theta\frac{\partial\left(\Phi\tilde{\kappa}\tilde{c}\right)}{\partial\tilde{c}}\divg\delta\mathbf{v}\right)\,dv\,,\label{eq250} \end{equation} \]
where
\[ \begin{equation} \boldsymbol{\sigma}_{c}^{\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\mathbf{S}^{e}}{\partial\tilde{c}}\cdot\mathbf{F}^{T}\label{eq251} \end{equation} \]
represents the spatial tangent of the stress with respect to the effective concentration. The next term is
\[ \begin{equation} D\left(\mathbf{W}\cdot\Grad\delta\tilde{p}\right)\left[\Delta\tilde{c}\right]\,dV=\grad\delta\tilde{p}\cdot\mathbf{w}_{c}^{\prime}\,dv\,,\label{eq252} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{c}^{\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{W}\left[\Delta\tilde{c}\right]=-\Delta\tilde{c}\,\tilde{\mathbf{k}}_{c}^{\prime}\cdot\left(\grad\tilde{p}+R\theta\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad\tilde{c}\right)\\ & -R\theta\tilde{\mathbf{k}}\cdot\left[\Delta\tilde{c}\left(\frac{\partial}{\partial\tilde{c}}\left(\frac{\tilde{\kappa}}{d_{0}}\right)\mathbf{d}+\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}_{c}^{\prime}\right)\cdot\grad\tilde{c}+\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad\Delta\tilde{c}\right]\,, \end{aligned} \label{eq253} \end{equation} \]
and
\[ \begin{equation} \tilde{\mathbf{k}}_{c}^{\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\tilde{\mathbf{K}}}{\partial\tilde{c}}\cdot\mathbf{F}^{T}\label{eq254} \end{equation} \]
is the spatial tangent of the effective hydraulic permeability with respect to the effective concentration.
The next term reduces to
\[ \begin{equation} -D\left(\delta\tilde{p}\frac{\partial J}{\partial t}\right)\left[\Delta\tilde{c}\right]\,dV=0\,.\label{eq255} \end{equation} \]
The following term is
\[ \begin{equation} D\left(\mathbf{J}\cdot\Grad\delta\tilde{c}\right)\left[\Delta\tilde{c}\right]\,dV=\grad\delta\tilde{c}\cdot\mathbf{j}_{c}^{\prime}\,dv\,,\label{eq256} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{j}_{c}^{\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{J}\left[\Delta\tilde{c}\right]\\ & =\Delta\tilde{c}\left(\frac{\partial\tilde{\kappa}}{\partial\tilde{c}}\mathbf{d}+\tilde{\kappa}\mathbf{d}_{c}^{\prime}\right)\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}\right)\\ & -\varphi^{w}\tilde{\kappa}\mathbf{d}\cdot\grad\Delta\tilde{c}+\tilde{\kappa}\frac{\tilde{c}}{d_{0}}\mathbf{d}\cdot\mathbf{w}_{c}^{\prime}\,, \end{aligned} \label{eq257} \end{equation} \]
and
\[ \begin{equation} \mathbf{d}_{c}^{\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\mathbf{D}}{\partial\tilde{c}}\cdot\mathbf{F}^{T}\label{eq258} \end{equation} \]
is the spatial tangent of the diffusivity with respect to the effective concentration.
The last term is
\[ \begin{equation} D\left(\frac{\partial\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)}{\partial t}\delta\tilde{c}\right)\left[\Delta\tilde{c}\right]\,dV=\delta\tilde{c}\frac{\varphi^{w}}{\Delta t}\frac{\partial\left(\tilde{\kappa}\tilde{c}\right)}{\partial\tilde{c}}\Delta\tilde{c}\,dv\,,\label{eq259} \end{equation} \]
where we similarly used a backward difference scheme to discretize the time derivative.
Linearization of External Virtual Work
The linearization of \(\delta W_{\text{ext}}\) in \eqref{eq226} depends on whether natural boundary conditions are prescribed as area densities or total net values over an area. Thus, in the case when \(\mathbf{t}\,da\) (net force), \(w_{n}da\) (net volumetric flow rate), or \(j_{n}da\) (net molar flow rate) are prescribed over the elemental area \(da\), there is no variation in \(\delta W_{\text{ext}}\) and it follows that \(D\delta W_{\text{ext}}=0\). Alternatively, in the case when \(\mathbf{t}\), \(w_{n}\) or \(j_{n}\) are prescribed, the linearization may be performed by evaluating the integral in the parametric space of the boundary surface \(\partial b\), with parametric coordinates \(\left(\eta^{1},\eta^{2}\right)\). Accordingly, for a point \(\mathbf{x}\left(\eta^{1},\eta^{2}\right)\) on \(\partial b\), surface tangents (covariant basis vectors) are given by
\[ \begin{equation} \mathbf{g}_{\alpha}=\frac{\partial\mathbf{x}}{\partial\eta^{\alpha}},\quad\left(\alpha=1,2\right)\label{eq260} \end{equation} \]
and the outward unit normal is
\[ \begin{equation} \mathbf{n}=\frac{\mathbf{g}_{1}\times\mathbf{g}_{2}}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,.\label{eq261} \end{equation} \]
The elemental area on \(\partial b\) is \(da=\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|d\eta^{1}d\eta^{2}\). Consequently, the external virtual work integral may be rewritten as
\[ \begin{equation} \delta W_{\text{ext}}=\int_{\partial b}\left(\delta\mathbf{v}\cdot\mathbf{t}+\delta\tilde{p}\,w_{n}+\delta\tilde{c}\,j_{n}\right)\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|d\eta^{1}d\eta^{2}\,.\label{eq262} \end{equation} \]
The directional derivative of \(\delta W_{\text{ext}}\) may then be applied directly to its integrand, since the parametric space is invariant .
If we restrict traction boundary conditions to the special case of normal tractions, then \(\mathbf{t}=t_{n}\mathbf{n}\) where \(t_{n}\) is the prescribed normal traction component. Then it can be shown that the linearization of \(\delta W_{\text{ext}}\) along \(\Delta\mathbf{u}\) produces
\[ \begin{equation} D\left(\delta W_{\text{ext}}\right)\left[\Delta\mathbf{u}\right]=\int_{\partial b}\left(t_{n}\delta\mathbf{v}+w_{n}\delta\tilde{p}\,\mathbf{n}+j_{n}\delta\tilde{c}\,\mathbf{n}\right)\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)d\eta^{1}d\eta^{2}\,.\label{eq263} \end{equation} \]
The linearizations along \(\Delta\tilde{p}\) and \(\Delta\tilde{c}\) reduce to zero, \(D\left(\delta W_{\text{ext}}\right)\left[\Delta\tilde{p}\right]=0\) and \(D\left(\delta W_{\text{ext}}\right)\left[\Delta\tilde{c}\right]=0\).
Discretization
To discretize the virtual work relations, let
\[ \begin{equation} \begin{aligned}\delta\mathbf{v} & =\sum\limits_{a=1}^{m}N_{a}\delta\mathbf{v}_{a}\,, & \Delta\mathbf{u} & =\sum\limits_{b=1}^{m}N_{b}\Delta\mathbf{u}_{b}\,,\\ \delta\tilde{p} & =\sum\limits_{a=1}^{m}N_{a}\delta\tilde{p}_{a}\,, & \Delta\tilde{p} & =\sum\limits_{b=1}^{m}N_{b}\Delta\tilde{p}_{b}\,,\\ \delta\tilde{c} & =\sum\limits_{a=1}^{m}N_{a}\delta\tilde{c}_{a}\,, & \Delta\tilde{c} & =\sum\limits_{b=1}^{m}N_{b}\Delta\tilde{c}_{b}\,, \end{aligned} \label{eq264} \end{equation} \]
where \(N_{a}\) represents the interpolation functions over an element, \(\delta\mathbf{v}_{a}\), \(\delta\tilde{p}_{a}\), \(\delta\tilde{c}_{a}\), \(\Delta\mathbf{u}_{a}\), \(\Delta\tilde{p}_{a}\) and \(\Delta\tilde{c}_{a}\) respectively represent the nodal values of \(\delta\mathbf{v}\), \(\delta\tilde{p}\), \(\delta\tilde{c}\), \(\Delta\mathbf{u}\), \(\Delta\tilde{p}\) and \(\Delta\tilde{c}\); \(m\) is the number of nodes in an element.
The discretized form of \(\delta W_{\text{int}}\) in \eqref{eq226} 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}{ccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{r}_{a}^{u}\\ r_{a}^{p}\\ r_{a}^{c} \end{array}\right]\,,\label{eq265} \end{equation} \]
where \(n_{e}\) is the number of elements in \(b\), \(n_{\text{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 current spatial configuration to the parametric space of the element. In the above expression,
\[ \begin{equation} \begin{aligned}\mathbf{r}_{a}^{u} & =\boldsymbol{\sigma}\cdot\grad N_{a}\,,\\ r_{a}^{p} & =\mathbf{w}\cdot\grad N_{a}-N_{a}\frac{1}{J}\frac{\partial J}{\partial t}\,,\\ r_{a}^{c} & =\mathbf{j}\cdot\grad N_{a}-N_{a}\frac{1}{J}\frac{\partial}{\partial t}\left(J\varphi^{w}\tilde{\kappa}\tilde{c}\right)\,, \end{aligned} \label{eq266} \end{equation} \]
and it is understood that \(J_{\eta}\), \(\mathbf{r}_{a}^{u}\), \(r_{a}^{p}\) and \(r_{a}^{c}\) are evaluated at the parametric coordinates of the \(k-\)th integration point. Since the parametric space is invariant, time derivatives are evaluated in a material frame. For example, the time derivative \(D^{s}J\left(\mathbf{x},t\right)/Dt\) appearing in \eqref{eq226} becomes \(\partial J\left(\eta_{k},t\right)/\partial t\) when evaluated at the parametric coordinates \(\eta_{k}=\left(\eta_{k}^{1},\eta_{k}^{2},\eta_{k}^{3}\right)\) of the \(k-\)th integration point.
Similarly, the discretized form of \(D\delta W_{\text{int}}=D\delta W_{\text{int}}\left[\Delta\mathbf{u}\right]+D\delta W_{\text{int}}\left[\Delta\tilde{p}\right]+D\delta W_{\text{int}}\left[\Delta\tilde{c}\right]\) 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}\sum\limits_{b=1}^{m}\left[\begin{array}{ccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}\end{array}\right]\cdot\left[\begin{array}{ccc} \mathbf{K}_{ab}^{uu} & \mathbf{k}_{ab}^{up} & \mathbf{k}_{ab}^{uc}\\ \mathbf{k}_{ab}^{pu} & k_{ab}^{pp} & k_{ab}^{pc}\\ \mathbf{k}_{ab}^{cu} & k_{ab}^{cp} & k_{ab}^{cc} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta\tilde{p}_{b}\\ \Delta\tilde{c}_{b} \end{array}\right]\,,\label{eq267} \end{equation} \]
where the terms in the first column are the discretized form of the linearization along \(\Delta\mathbf{u}\):
\[ \begin{equation} \mathbf{K}_{ab}^{uu}=\grad N_{a}\cdot\boldsymbol{\mathcal{C}}\cdot\grad N_{b}+\left(\grad N_{a}\cdot\boldsymbol{\sigma}\cdot\grad N_{b}\right)\mathbf{I}\,,\label{eq268} \end{equation} \]
\[ \begin{equation} \mathbf{k}_{ab}^{pu}=\left(\mathbf{w}_{b}^{u}\right)^{T}\cdot\grad N_{a}+N_{a}\mathbf{q}_{b}^{pu}\,,\label{eq269} \end{equation} \]
\[ \begin{equation} \mathbf{k}_{ab}^{cu}=\left(\mathbf{j}_{b}^{u}\right)^{T}\cdot\grad N_{a}+N_{a}\mathbf{q}_{b}^{cu}\,,\label{eq270} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{j}_{b}^{u} & =J\frac{\partial\tilde{\kappa}}{\partial J}\left[\mathbf{d}\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}\right)\right]\otimes\grad N_{b}+\tilde{\kappa}\left({-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}}\right)\cdot\mathbf{d}\cdot\grad N_{b}\\ & +\tilde{\kappa}\left({-\varphi^{s}\left({\mathbf{d}\cdot\grad\tilde{c}}\right)\otimes\grad N_{b}+\frac{\tilde{c}}{d_{0}}\left[{2\left({\grad N_{b}\cdot\mathbf{w}}\right)\mathbf{d}-\left({\mathbf{d}\cdot\mathbf{w}}\right)\otimes\grad N_{b}}\right]}\right)+\tilde{\kappa}\frac{\tilde{c}}{d_{0}}\mathbf{d}\cdot\mathbf{w}_{b}^{u}\;, \end{aligned} \label{eq272} \end{equation} \]
\[ \begin{equation} \mathbf{q}_{b}^{pu}=-\frac{1}{\Delta t}\grad N_{b},\label{eq273} \end{equation} \]
\[ \begin{equation} \mathbf{q}_{b}^{cu}=\tilde{c}\frac{\partial\left({J\phi^{w}\tilde{\kappa}}\right)}{\partial J}\mathbf{q}_{b}^{pu}.\label{eq274} \end{equation} \]
The terms in the second column of the stiffness matrix in \eqref{eq267} are the discretized form of the linearization along \(\Delta\tilde{p}\):
\[ \begin{equation} \mathbf{k}_{ab}^{up}=-N_{b}\grad N_{a},\label{eq275} \end{equation} \]
\[ \begin{equation} k_{ab}^{pp}=-\grad N_{a}\cdot\tilde{\mathbf{k}}\cdot\grad N_{b},\label{eq276} \end{equation} \]
\[ \begin{equation} k_{ab}^{cp}=-\frac{\tilde{\kappa}\tilde{c}}{d_{0}}\grad N_{a}\cdot\mathbf{d}\cdot\tilde{\mathbf{k}}\cdot\grad N_{b}.\label{eq277} \end{equation} \]
The terms in the third column of the stiffness matrix in \eqref{eq267} are the discretized form of the linearization along \(\Delta\tilde{c}\):
\[ \begin{equation} \mathbf{k}_{ab}^{uc}=N_{b}\left(\boldsymbol{\sigma}_{c}^{\prime}\cdot\grad N_{a}-R\theta\frac{\partial\left(\Phi\tilde{\kappa}\tilde{c}\right)}{\partial\tilde{c}}\grad N_{a}\right)\,,\label{eq278} \end{equation} \]
\[ \begin{equation} k_{ab}^{pc}=\grad N_{a}\cdot\mathbf{w}_{b}^{c}\,,\label{eq279} \end{equation} \]
\[ \begin{equation} k_{ab}^{cc}=\grad N_{a}\cdot\mathbf{j}_{b}^{c}+N_{a}q_{b}^{c}\,,\label{eq280} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{b}^{c} & =-N_{b}\tilde{\mathbf{k}}_{c}^{\prime}\cdot\left(\grad\tilde{p}+R\theta\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad\tilde{c}\right)\\ & -R\theta\tilde{\mathbf{k}}\cdot\left[N_{b}\left(\frac{\partial}{\partial\tilde{c}}\left(\frac{\tilde{\kappa}}{d_{0}}\right)\mathbf{d}+\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}_{c}^{\prime}\right)\cdot\grad\tilde{c}+\frac{\tilde{\kappa}}{d_{0}}\mathbf{d}\cdot\grad N_{b}\right]\,, \end{aligned} \label{eq281} \end{equation} \]
\[ \begin{equation} \mathbf{j}_{b}^{c}=N_{b}\left(\frac{\partial\tilde{\kappa}}{\partial\tilde{c}}\mathbf{d}+\tilde{\kappa}\mathbf{d}_{c}^{\prime}\right)\cdot\left(-\varphi^{w}\grad\tilde{c}+\frac{\tilde{c}}{d_{0}}\mathbf{w}\right)+\tilde{\kappa}\mathbf{d}\cdot\left(-\varphi^{w}\grad N_{b}+\frac{\tilde{c}}{d_{0}}\mathbf{w}_{b}^{c}\right)\,,\label{eq282} \end{equation} \]
\[ \begin{equation} q_{b}^{c}=-N_{b}\frac{\phi^{w}}{\Delta t}\frac{\partial\left(\tilde{\kappa}\tilde{c}\right)}{\partial\tilde{c}}\,.\label{eq283} \end{equation} \]
The discretization of \(\delta W_{\text{ext}}\) in \eqref{eq226} has the form
\[ \begin{equation} \delta W_{\text{ext}}=\sum\limits_{e=1}^{n_{e}}\sum\limits_{k=1}^{n_{\mbox{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\left[\begin{array}{ccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}\end{array}\right]\cdot\left[\begin{array}{c} N_{a}t_{n}\mathbf{n}\\ N_{a}w_{n}\\ N_{a}j_{n} \end{array}\right]\,,\label{eq284} \end{equation} \]
where \(J_{\eta}=\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\). The summation is performed over all surface elements on which these boundary conditions are prescribed. The discretization of \(-D\delta W_{\text{ext}}\) has the form
\[ \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}\sum\limits_{b=1}^{m}\left[\begin{array}{ccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}\end{array}\right]\cdot\left[\begin{array}{ccc} \mathbf{K}_{ab}^{uu} & \mathbf{0} & \mathbf{0}\\ \mathbf{k}_{ab}^{pu} & 0 & 0\\ \mathbf{k}_{ab}^{cu} & 0 & 0 \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta\tilde{p}_{b}\\ \Delta\tilde{c}_{b} \end{array}\right]\,,\label{eq285} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{K}_{ab}^{uu} & =t_{n}N_{a}\boldsymbol{\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}^{pu} & =-w_{n}N_{a}\left(\frac{\partial N_{b}}{\partial\eta^{1}}\mathbf{g}_{2}-\frac{\partial N_{b}}{\partial\eta^{2}}\mathbf{g}_{1}\right)\times\mathbf{n}\,,\\ \mathbf{k}_{ab}^{cu} & =-j_{n}N_{a}\left(\frac{\partial N_{b}}{\partial\eta^{1}}\mathbf{g}_{2}-\frac{\partial N_{b}}{\partial\eta^{2}}\mathbf{g}_{1}\right)\times\mathbf{n}\,. \end{aligned} \label{eq286} \end{equation} \]
In this expression, \(\boldsymbol{\mathcal{A}}\left\{ \mathbf{v}\right\}\) is the antisymmetric tensor whose dual vector is \(\mathbf{v}\) (such that \(\boldsymbol{\mathcal{A}}\left\{ \mathbf{v}\right\} \cdot\mathbf{q}=\mathbf{v}\times\mathbf{q}\) for any vector \(\mathbf{q})\).