The virtual work integral for a mixture of intrinsically incompressible constituents combines the balance of momentum for the mixture, the balance of mass for the mixture, and the balance of mass for each of the solutes. In addition, for charged mixtures, the condition of (2.9-8) may be enforced as a penalty constraint on each solute mass balance equation:
\[ \begin{equation} \begin{aligned}\delta W & =\int_{b}\delta\mathbf{v}\cdot\mbox{div}\boldsymbol{\sigma}\,dv\\ & +\int_{b}\delta\tilde{p}\,\mbox{div}\left(\mathbf{v}^{s}+\mathbf{w}\right)\,dv\\ & +\sum\limits_{\alpha\ne s,w}\int_{b}\delta\tilde{c}^{\alpha}\left[\frac{1}{J^{s}}\frac{D^{s}}{Dt}\left(J^{s}\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)+\mbox{div}\mathbf{j}^{\alpha}+\sum\limits_{\beta\ne s,w}z^{\beta}\mbox{div}\mathbf{j}^{\beta}\right]\,dv\,, \end{aligned} \label{eq287} \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}^{\alpha}\) is the virtual molar energy of solute \(\alpha\). Here, \(b\) represents the mixture domain in the spatial frame and \(dv\) is an elemental volume in \(b\). Applying the divergence theorem, \(\delta W\) may be split into internal and external contributions to the virtual work, \(\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}\,dv+\int_{b}\left(\mathbf{w}\cdot\mbox{grad}\delta\tilde{p}-\frac{\delta\tilde{p}}{J^{s}}\frac{D^{s}J^{s}}{Dt}\right)\,dv\\ & +\sum\limits_{\alpha\ne s,w}\int_{b}\left[\mathbf{j}^{\alpha}\cdot\mbox{grad}\delta\tilde{c}^{\alpha}-\frac{\delta\tilde{c}^{\alpha}}{J^{s}}\frac{D^{s}}{Dt}\left(J^{s}\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)\right]\,dv\\ & +\sum\limits_{\alpha\ne s,w}\int_{b}\mbox{grad}\delta\tilde{c}^{\alpha}\cdot\sum\limits_{\beta\ne s,w}z^{\beta}\mathbf{j}^{\beta}\,dv\,, \end{aligned} \label{eq288} \end{equation} \]
and
\[ \begin{equation} \delta W_{\text{ext}}=\int_{\partial b}\left[\delta\mathbf{v}\cdot\mathbf{t}+\delta\tilde{p}\,w_{n}+\sum\limits_{\alpha\ne s,w}\delta\tilde{c}^{\alpha}\left(j_{n}^{\alpha}+\sum\limits_{\beta\ne s,w}z^{\beta}j_{n}^{\beta}\right)\right]\,da\,.\label{eq289} \end{equation} \]
In these expressions, \(\delta\mathbf{D}=\left(\mbox{grad}\delta\mathbf{v}+\mbox{grad}^{T}\delta\mathbf{v}\right)/2\), \(\partial b\) is the boundary of \(b\), and \(da\) is an elemental area on \(\partial b\). In this finite element formulation, \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}^{\alpha}\) are used as nodal variables, and essential boundary conditions may be prescribed on these variables. Natural boundary conditions are prescribed to the mixture traction, \(\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n}\), normal fluid flux, \(w_{n}=\mathbf{w}\cdot\mathbf{n}\), and normal solute flux, \(j_{n}^{\alpha}=\mathbf{j}^{\alpha}\cdot\mathbf{n}\), where \(\mathbf{n}\) is the outward unit normal to \(\partial b\). To solve the system \(\delta W=0\) for nodal values of \(\mathbf{u}\), \(\tilde{p}\) and \(\tilde{c}^{\alpha}\), it is necessary to linearize these equations, as shown for example in Sections Linearization of Internal Virtual Work-Linearization of External Virtual Work for biphasic-solute materials. If the mixture is charged, it is also necessary to solve for the electric potential \(\psi\) by solving the algebraic relation of the electroneutrality condition in (2.9-4), which may be rewritten as
\[ \begin{equation} c^{F}+\sum\limits_{\beta\ne s,w}z^{\beta}\tilde{\kappa}^{\beta}\tilde{c}^{\beta}=0\,.\label{eq290} \end{equation} \]
In the special case of a triphasic mixture, where solutes consist of two counter-ions (\(\alpha=+,-)\), this equation may be solved in closed form to produce
\[ \begin{equation} \psi=\frac{1}{z^{\alpha}}\frac{R\theta}{F_{c}}\ln\left(\frac{2z^{\alpha}\hat{\kappa}^{\alpha}\tilde{c}^{\alpha}}{-c^{F}\pm\sqrt{\left(c^{F}\right)^{2}+4\left(z^{\alpha}\right)^{2}\left(\hat{\kappa}^{+}\tilde{c}^{+}\right)\left(\hat{\kappa}^{-}\tilde{c}^{-}\right)}}\right),\quad\alpha=+,-,\label{eq291} \end{equation} \]
Only the positive root is valid in the argument of the logarithm function.
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\mathbf{\dot{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{eq292} \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\sum\limits_{\beta}\tilde{\kappa}^{\beta}\tilde{c}^{\beta}\right)\left(\mathbf{I}\otimes\mathbf{I}-2\mathbf{I}\odot\mathbf{I}\right)-R\theta\sum\limits_{\beta}\tilde{c}^{\beta}\,J\frac{\partial\left(\Phi\tilde{\kappa}^{\beta}\right)}{\partial J}\mathbf{I}\otimes\mathbf{I}\,,\label{eq293} \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}\oslash\mathbf{F}\right):2\frac{\partial\mathbf{S}^{e}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right)\,.\label{eq294} \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{eq295} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{u}^{\prime} & \equiv-\left(\tilde{\boldsymbol{\mathcal{K}}}:\Delta\boldsymbol{\varepsilon}\right)\cdot\left(\grad\tilde{p}+R\theta\sum\limits_{\beta}\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\mathbf{d}^{\beta}\cdot\grad\tilde{c}^{\beta}\right)\\ & -\tilde{\mathbf{k}}\cdot R\theta\sum\limits_{\beta}\left(\left[J\frac{\partial}{\partial J}\left(\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\right)-\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\right]\left(\divg\Delta\mathbf{u}\right)\mathbf{d}^{\beta}+\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\left(2\Delta\boldsymbol{\varepsilon}\cdot\mathbf{d}^{\beta}+\boldsymbol{\mathcal{D}}^{\beta}:\Delta\boldsymbol{\varepsilon}\right)\right)\cdot\grad\tilde{c}^{\beta} \end{aligned} \label{eq296} \end{equation} \]
with
\[ \begin{equation} \begin{aligned}\tilde{\boldsymbol{\mathcal{K}}} & =J^{-1}\left(\mathbf{F}\oslash\mathbf{F}\right):2\frac{\partial\tilde{\mathbf{K}}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right)\,,\\ \boldsymbol{\mathcal{D}}^{\alpha} & =J^{-1}\left(\mathbf{F}\oslash\mathbf{F}\right):2\frac{\partial\mathbf{D}^{\alpha}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right)\,, \end{aligned} \label{eq297} \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 (3.3-8)\(_{1}\), the evaluation of \(\tilde{\boldsymbol{\mathcal{K}}}\) is rather involved and it can be shown that
\[ \begin{equation} \tilde{\boldsymbol{\mathcal{K}}}=\left(\tilde{\mathbf{k}}\oslash\tilde{\mathbf{k}}\right):\left[\left(\mathbf{k}^{-1}\oslash\mathbf{k}^{-1}\right):\boldsymbol{\mathcal{K}}+\boldsymbol{\mathcal{G}}\right]\,,\label{eq298} \end{equation} \]
where
\[ \begin{equation} \boldsymbol{\mathcal{K}}=J^{-1}\left(\mathbf{F}\oslash\mathbf{F}\right):2\frac{\partial\mathbf{K}}{\partial\mathbf{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right)\,,\label{eq300} \end{equation} \]
and
\[ \begin{equation} \boldsymbol{\mathcal{G}}=\frac{R\theta}{\varphi^{w}}\sum_{\alpha}\frac{\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\left(\begin{aligned}\left(\frac{1}{\varphi^{w}}-\frac{J}{\tilde{\kappa}^{\alpha}}\frac{\partial\tilde{\kappa}^{\alpha}}{\partial J}\right)\left(\mathbf{I}-\frac{\mathbf{d}^{\alpha}}{d_{0}^{\alpha}}\right)\otimes\mathbf{I}\\ +2\left(\mathbf{I}\odot\frac{\mathbf{d}^{\alpha}}{d_{0}^{\alpha}}+\frac{\mathbf{d}^{\alpha}}{d_{0}^{\alpha}}\odot\mathbf{I}-\mathbf{I}\odot\mathbf{I}\right)\\ -\frac{\mathbf{d}^{\alpha}}{d_{0}^{\alpha}}\otimes\mathbf{I}+\frac{\boldsymbol{\mathcal{D}}^{\alpha}}{d_{0}^{\alpha}} \end{aligned} \right)\,.\label{eq299} \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{eq301} \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{eq302} \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}^{\alpha}\cdot\Grad\delta\tilde{c}^{\alpha}\right)\left[\Delta\mathbf{u}\right]\,dV=\grad\delta\tilde{c}^{\alpha}\cdot\mathbf{j}_{u}^{\alpha\prime}\,dv\,,\label{eq303} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{j}_{u}^{\alpha\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{J}\left[\Delta\mathbf{u}\right]=\left(J\frac{\partial\tilde{\kappa}^{\alpha}}{\partial J}\left(\divg\Delta\mathbf{u}\right)\mathbf{d}^{\alpha}+\tilde{\kappa}^{\alpha}\boldsymbol{\mathcal{D}}^{\alpha}:\Delta\boldsymbol{\varepsilon}\right)\cdot\mathbf{g}^{\alpha}\\ & +\tilde{\kappa}^{\alpha}\mathbf{d}^{\alpha}\cdot\left(-\varphi^{s}\left(\divg\Delta\mathbf{u}\right)\grad\tilde{c}^{\alpha}+\left(2\Delta\boldsymbol{\varepsilon}-\left(\divg\Delta\mathbf{u}\right)\mathbf{I}\right)\cdot\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{w}+\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{w}_{u}^{\prime}\right) \end{aligned} \,.\label{eq304} \end{equation} \]
where
\[ \begin{equation} \mathbf{g}^{\alpha}=-\varphi^{w}\grad\tilde{c}^{\alpha}+\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{w}\,.\label{eq305} \end{equation} \]
Using a backward difference scheme for the time derivative, the last term is
\[ \begin{equation} D\left(\delta\tilde{c}^{\alpha}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)}{\partial t}\right)\left[\Delta\mathbf{u}\right]\,dV=\frac{\delta\tilde{c}^{\alpha}}{\Delta t}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\right)}{\partial J}\tilde{c}^{\alpha}\divg\Delta\mathbf{u}\,dv\,.\label{eq306} \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{eq307} \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{eq308} \end{equation} \]
\[ \begin{equation} D\left(\mathbf{J}^{\alpha}\cdot\Grad\delta\tilde{c}^{\alpha}-\delta\tilde{c}^{\alpha}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)}{\partial t}\right)\left[\Delta\tilde{p}\right]\,dV=-\frac{\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\grad\delta\tilde{c}^{\alpha}\cdot\mathbf{d}^{\alpha}\cdot\tilde{\mathbf{k}}\cdot\grad\Delta\tilde{p}\,dv\,.\label{eq309} \end{equation} \]
Linearization along \(\Delta\tilde{c}^{\gamma}\)
The linearization of the first term in \(\delta W_{\text{int}}\) along \(\Delta\tilde{c}^{\gamma}\) yields
\[ \begin{equation} D\left(\mathbf{S}:\delta\mathbf{\dot{E}}\right)\left[\Delta\tilde{c}^{\gamma}\right]dV=\Delta\tilde{c}^{\gamma}\left(\boldsymbol{\sigma}_{\gamma}^{\prime}:\delta\mathbf{d}-R\theta\frac{\partial\left(\Phi\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)}{\partial\tilde{c}^{\gamma}}\divg\delta\mathbf{v}\right)\,dv\,,\label{eq310} \end{equation} \]
where
\[ \begin{equation} \boldsymbol{\sigma}_{\gamma}^{\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\mathbf{S}^{e}}{\partial\tilde{c}^{\gamma}}\cdot\mathbf{F}^{T}\label{eq311} \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}^{\gamma}\right]\,dV=\grad\delta\tilde{p}\cdot\mathbf{w}_{\gamma}^{\prime}\,dv\,,\label{eq312} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{\gamma}^{\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{W}\left[\Delta\tilde{c}^{\gamma}\right]=-\tilde{\mathbf{k}}_{\gamma}^{\prime}\cdot\left(\grad\tilde{p}+R\theta\sum\limits_{\beta}\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\mathbf{d}^{\beta}\cdot\grad\tilde{c}^{\beta}\right)\Delta\tilde{c}^{\gamma}\\ & -R\theta\tilde{\mathbf{k}}\cdot\frac{\tilde{\kappa}^{\gamma}}{d_{0}^{\gamma}}\mathbf{d}^{\gamma}\cdot\grad\Delta\tilde{c}^{\gamma}-\Delta\tilde{c}^{\gamma}R\theta\tilde{\mathbf{k}}\cdot\sum\limits_{\beta}\left(\frac{\partial}{\partial\tilde{c}^{\gamma}}\left(\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\right)\mathbf{d}^{\beta}+\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\mathbf{d}_{\gamma}^{\beta\prime}\right)\cdot\grad\tilde{c}^{\beta} \end{aligned} \label{eq313} \end{equation} \]
and
\[ \begin{equation} \tilde{\mathbf{k}}_{\gamma}^{\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\tilde{\mathbf{K}}}{\partial\tilde{c}^{\gamma}}\cdot\mathbf{F}^{T},\quad\mathbf{d}_{\gamma}^{\beta\prime}=J^{-1}\mathbf{F}\cdot\frac{\partial\mathbf{d}^{\beta}}{\partial\tilde{c}^{\gamma}}\cdot\mathbf{F}^{T}\label{eq314} \end{equation} \]
are the spatial tangents of the effective hydraulic permeability and solute diffusivity 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}^{\gamma}\right]\,dV=0\,.\label{eq315} \end{equation} \]
The following term is
\[ \begin{equation} D\left(\mathbf{J}^{\alpha}\cdot\Grad\delta\tilde{c}^{\alpha}\right)\left[\Delta\tilde{c}^{\gamma}\right]\,dV=\grad\delta\tilde{c}^{\alpha}\cdot\mathbf{j}_{\gamma}^{\alpha\prime}\,dv\,,\label{eq316} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{j}_{\gamma}^{\alpha\prime} & \equiv J^{-1}\mathbf{F}\cdot D\mathbf{J}^{\alpha}\left[\Delta\tilde{c}^{\gamma}\right]=\Delta\tilde{c}^{\gamma}\left(\frac{\partial\tilde{\kappa}^{\alpha}}{\partial\tilde{c}^{\gamma}}\mathbf{d}^{\alpha}+\tilde{\kappa}^{\alpha}\mathbf{d}_{\gamma}^{\alpha\prime}\right)\cdot\mathbf{g}^{\alpha}\\ & +\tilde{\kappa}^{\alpha}\mathbf{d}^{\alpha}\cdot\left(-\varphi^{w}\delta_{\alpha\gamma}\grad\Delta\tilde{c}^{\gamma}+\frac{\Delta\tilde{c}^{\gamma}}{d_{0}^{\alpha}}\left(\delta_{\alpha\gamma}-\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\frac{\partial d_{0}^{\alpha}}{\partial\tilde{c}^{\gamma}}\right)\mathbf{w}+\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{w}_{\gamma}^{\prime}\right)\,. \end{aligned} \label{eq317} \end{equation} \]
The last term is
\[ \begin{equation} D\left(\frac{\partial\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)}{\partial t}\delta\tilde{c}^{\alpha}\right)\left[\Delta\tilde{c}^{\gamma}\right]\,dV=\delta\tilde{c}^{\alpha}\frac{\varphi^{w}}{\Delta t}\frac{\partial\left(\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)}{\partial\tilde{c}^{\gamma}}\Delta\tilde{c}^{\gamma}\,dv\,,\label{eq318} \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{eq289} 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 \(\tilde{j}_{n}^{\alpha}da\) (net effective 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 \(\tilde{j}_{n}^{\alpha}\) 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{eq319} \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{eq320} \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}+\sum\limits_{\alpha\ne s,w}\delta\tilde{c}^{\alpha}\tilde{j}_{n}^{\alpha}\right)\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\,d\eta^{1}d\eta^{2}\,,\label{eq321} \end{equation} \]
where
\[ \begin{equation} \tilde{j}_{n}^{\alpha}=j_{n}^{\alpha}+\sum\limits_{\beta\ne s,w}z^{\beta}j_{n}^{\beta}\,.\label{eq322} \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}+\sum\limits_{\alpha\ne s,w}\delta\tilde{c}^{\alpha}\tilde{j}_{n}^{\alpha}\,\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{eq323} \end{equation} \]
The linearizations along \(\Delta\tilde{p}\) and \(\Delta\tilde{c}^{\gamma}\) 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}^{\gamma}\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}^{\alpha} & =\sum\limits_{a=1}^{m}N_{a}\delta\tilde{c}_{a}^{\alpha}\,, & \Delta\tilde{c}^{\gamma} & =\sum\limits_{b=1}^{m}N_{b}\Delta\tilde{c}_{b}^{\gamma}\,, \end{aligned} \label{eq324} \end{equation} \]
where \(N_{a}\) represents the interpolation functions over an element, \(\delta\mathbf{v}_{a}\), \(\delta\tilde{p}_{a}\), \(\delta\tilde{c}_{a}^{\alpha}\), \(\Delta\mathbf{u}_{a}\), \(\Delta\tilde{p}_{a}\) and \(\Delta\tilde{c}_{a}^{\gamma}\) respectively represent the nodal values of \(\delta\mathbf{v}\), \(\delta\tilde{p}\), \(\delta\tilde{c}^{\alpha}\), \(\Delta\mathbf{u}\), \(\Delta\tilde{p}\) and \(\Delta\tilde{c}^{\gamma}\); \(m\) is the number of nodes in an element.
The discretized form of \(\delta W_{\text{int}}\) in (3.3-3) 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}{cccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}^{\alpha} & \delta\tilde{c}_{a}^{\beta}\end{array}\right]\cdot\left[\begin{array}{c} \mathbf{r}_{a}^{u}\\ r_{a}^{p}\\ r_{a}^{\alpha}\\ r_{a}^{\beta} \end{array}\right]\,,\label{eq325} \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}^{\alpha} & =\mathbf{j}^{\alpha}\cdot\grad N_{a}-N_{a}\frac{1}{J}\frac{\partial}{\partial t}\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}\right)\,, \end{aligned} \label{eq326} \end{equation} \]
and it is understood that \(J_{\eta}\), \(\mathbf{r}_{a}^{u}\), \(r_{a}^{p}\) and \(r_{a}^{\alpha}\) 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 (3.3-3) 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. All time derivatives are discretized using a backward difference scheme.
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]+\sum\nolimits_{\gamma}D\delta W_{\text{int}}\left[\Delta\tilde{c}^{\gamma}\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}{cccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}^{\alpha} & \delta\tilde{c}_{a}^{\beta}\end{array}\right]\cdot\left[\begin{array}{cccc} \mathbf{K}_{ab}^{uu} & \mathbf{k}_{ab}^{up} & \mathbf{k}_{ab}^{u\alpha} & \mathbf{k}_{ab}^{u\beta}\\ \mathbf{k}_{ab}^{pu} & k_{ab}^{pp} & k_{ab}^{p\alpha} & k_{ab}^{p\beta}\\ \mathbf{k}_{ab}^{\alpha u} & k_{ab}^{\alpha p} & k_{ab}^{\alpha\alpha} & k_{ab}^{\alpha\beta}\\ \mathbf{k}_{ab}^{\beta u} & k_{ab}^{\beta p} & k_{ab}^{\beta\alpha} & k_{ab}^{\beta\beta} \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta\tilde{p}_{b}\\ \Delta\tilde{c}_{b}^{\alpha}\\ \Delta\tilde{c}_{b}^{\beta} \end{array}\right]\,,\label{eq327} \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{eq328} \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{eq329} \end{equation} \]
\[ \begin{equation} \mathbf{k}_{ab}^{\alpha u}=\left(\mathbf{j}_{b}^{\alpha u}+\sum\nolimits_{\beta}z^{\beta}\mathbf{j}_{b}^{\beta u}\right)^{T}\cdot\grad N_{a}+N_{a}\mathbf{q}_{b}^{\alpha u}\,,\label{eq330} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{b}^{u} & =\mathbf{g}^{p}\cdot\tilde{\boldsymbol{\mathcal{K}}}\cdot\grad N_{b}\\ & -R\theta\sum_{\beta}\frac{1}{d_{0}^{\beta}}\left(J\frac{\partial\tilde{\kappa}^{\beta}}{\partial J}-\tilde{\kappa}^{\beta}\right)\tilde{\mathbf{k}}\cdot\mathbf{d}^{\beta}\cdot\left(\grad\tilde{c}^{\beta}\otimes\grad N_{b}\right)\\ & -R\theta\sum_{\beta}\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\left[\tilde{\mathbf{k}}\cdot\left(\grad N_{b}\otimes\grad\tilde{c}^{\beta}\right)\cdot\mathbf{d}^{\beta}+\left(\grad N_{b}\cdot\mathbf{d}^{\beta}\cdot\grad\tilde{c}^{\beta}\right)\tilde{\mathbf{k}}\right]\\ & -R\theta\sum_{\beta}\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\tilde{\mathbf{k}}\cdot\left(\grad\tilde{c}^{\beta}\cdot\boldsymbol{\mathcal{D}}^{\beta}\cdot\grad N_{b}\right)\,, \end{aligned} \label{eq330b} \end{equation} \]
\[ \begin{equation} \mathbf{g}^{p}=-\grad\tilde{p}-R\theta\sum\limits_{\beta}\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\mathbf{d}^{\beta}\cdot\grad\tilde{c}^{\beta}\,.\label{eq344} \end{equation} \]
\[ \begin{equation} \begin{aligned}\mathbf{j}_{b}^{\alpha u} & =J\frac{\partial\tilde{\kappa}^{\alpha}}{\partial J}\mathbf{d}^{\alpha}\cdot\left(\mathbf{g}^{\alpha}\otimes\grad N_{b}\right)+\tilde{\kappa}^{\alpha}\mathbf{g}^{\alpha}\cdot\boldsymbol{\mathcal{D}}^{\alpha}\cdot\grad N_{b}\\ & +\tilde{\kappa}^{\alpha}\mathbf{d}^{\alpha}\cdot\left(\grad N_{b}\otimes\mathbf{w}-\mathbf{w}\otimes\grad N_{b}+\left(\grad N_{b}\cdot\mathbf{w}\right)\mathbf{I}+\mathbf{w}_{b}^{u}\right)\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\\ & -\varphi^{s}\tilde{\kappa}^{\alpha}\mathbf{d}^{\alpha}\cdot\grad\tilde{c}^{\alpha}\otimes\grad N_{b}\,, \end{aligned} \label{eq331} \end{equation} \]
\[ \begin{equation} \mathbf{g}^{\alpha}=-\varphi^{w}\grad\tilde{c}^{\alpha}+\frac{\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{w}\,,\label{eq331b} \end{equation} \]
\[ \begin{equation} \mathbf{q}_{b}^{pu}=-\frac{1}{\Delta t}\grad N_{b}\,,\label{eq332} \end{equation} \]
\[ \begin{equation} \mathbf{q}_{b}^{\alpha u}=\tilde{c}^{\alpha}\frac{\partial\left(J\varphi^{w}\tilde{\kappa}^{\alpha}\right)}{\partial J}\mathbf{q}_{b}^{pu}\,.\label{eq333} \end{equation} \]
The terms in the second column of the stiffness matrix in (3.3-44) are the discretized form of the linearization along \(\Delta\tilde{p}\):
\[ \begin{equation} \mathbf{k}_{ab}^{up}=-N_{b}\grad N_{a}\,,\label{eq334} \end{equation} \]
\[ \begin{equation} k_{ab}^{pp}=-\grad N_{a}\cdot\tilde{\mathbf{k}}\cdot\grad N_{b}\,,\label{eq335} \end{equation} \]
\[ \begin{equation} k_{ab}^{\alpha p}=\grad N_{a}\cdot\left(\mathbf{j}_{b}^{\alpha p}+\sum\limits_{\beta}z^{\beta}\mathbf{j}_{b}^{\beta p}\right)\,,\label{eq336} \end{equation} \]
where
\[ \begin{equation} \mathbf{j}_{b}^{\alpha p}=-\frac{\tilde{\kappa}^{\alpha}\tilde{c}^{\alpha}}{d_{0}^{\alpha}}\mathbf{d}^{\alpha}\cdot\tilde{\mathbf{k}}\cdot\grad N_{b}\,.\label{eq337} \end{equation} \]
The terms in the third column of the stiffness matrix in (3.3-44) are the discretized form of the linearization along \(\Delta\tilde{c}^{\gamma}\):
\[ \begin{equation} \mathbf{k}_{ab}^{u\alpha}=N_{b}\left(\boldsymbol{\sigma}_{\alpha}^{\prime}-R\theta\left[\Phi\tilde{\kappa}^{\alpha}+\sum\limits_{\beta}\left(\frac{\partial\Phi}{\partial\tilde{c}^{\alpha}}\tilde{\kappa}^{\beta}+\Phi\frac{\partial\tilde{\kappa}^{\beta}}{\partial\tilde{c}^{\alpha}}\right)\tilde{c}^{\beta}\right]\mathbf{I}\right)\cdot\grad N_{a}\,,\label{eq338} \end{equation} \]
\[ \begin{equation} k_{ab}^{p\alpha}=\grad N_{a}\cdot\mathbf{w}_{b}^{\alpha}\,,\label{eq339} \end{equation} \]
\[ \begin{equation} k_{ab}^{\alpha\gamma}=\grad N_{a}\cdot\left(\mathbf{j}_{b}^{\alpha\gamma}+\sum\limits_{\beta}z^{\beta}\mathbf{j}_{b}^{\beta\gamma}\right)+N_{a}q_{b}^{\alpha\gamma}\,,\label{eq340} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{w}_{b}^{\gamma} & =N_{b}\left(\tilde{\mathbf{k}}_{\gamma}^{\prime}\cdot\mathbf{g}^{p}-R\theta\tilde{\mathbf{k}}\cdot\sum\limits_{\beta}\left(\frac{\partial}{\partial\tilde{c}^{\gamma}}\left(\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\right)\mathbf{d}^{\beta}+\frac{\tilde{\kappa}^{\beta}}{d_{0}^{\beta}}\mathbf{d}_{c}^{\beta\gamma}\right)\cdot\grad\tilde{c}^{\beta}\right)\\ & -R\theta\tilde{\mathbf{k}}\cdot\frac{\tilde{\kappa}^{\gamma}}{d_{0}^{\gamma}}\mathbf{d}^{\gamma}\cdot\grad N_{b}\,, \end{aligned} \label{eq341} \end{equation} \]
and
\[ \begin{equation} \begin{aligned}\mathbf{j}_{b}^{\alpha\gamma} & =N_{b}\left(\frac{\partial\tilde{\kappa}^{\alpha}}{\partial\tilde{c}^{\gamma}}\mathbf{d}^{\alpha}+\tilde{\kappa}^{\alpha}\mathbf{d}_{c}^{\alpha\gamma}\right)\cdot\mathbf{g}^{\alpha}\\ & +\frac{\tilde{\kappa}^{\alpha}}{d_{0}^{\alpha}}\mathbf{d}^{\alpha}\cdot\left[\delta_{\alpha\gamma}\left(N_{b}\mathbf{w}-\varphi^{w}d_{0}^{\alpha}\grad N_{b}\right)+\tilde{c}^{\alpha}\left(\mathbf{w}_{b}^{\gamma}-\frac{N_{b}}{d_{0}^{\alpha}}\frac{\partial d_{0}^{\alpha}}{\partial\tilde{c}^{\gamma}}\mathbf{w}\right)\right]\,, \end{aligned} \label{eq342} \end{equation} \]
\[ \begin{equation} q_{b}^{\alpha\gamma}=-N_{b}\frac{\varphi^{w}}{\Delta t}\left(\frac{\partial\tilde{\kappa}^{\alpha}}{\partial\tilde{c}^{\gamma}}\tilde{c}^{\alpha}+\delta_{\alpha\gamma}\tilde{\kappa}^{\alpha}\right)\,.\label{eq343} \end{equation} \]
The discretization of \(\delta W_{\text{ext}}\) in \eqref{eq321} has the form
\[ \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}{cccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}^{\alpha} & \delta\tilde{c}_{a}^{\beta}\end{array}\right]\cdot\left[\begin{array}{c} N_{a}t_{n}\mathbf{n}\\ N_{a}w_{n}\\ N_{a}\tilde{j}_{n}^{\alpha}\\ N_{a}\tilde{j}_{n}^{\beta} \end{array}\right]\,,\label{eq345} \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_{\mbox{int}}^{\left(e\right)}}W_{k}J_{\eta}\sum\limits_{a=1}^{m}\sum\limits_{b=1}^{m}\left[\begin{array}{cccc} \delta\mathbf{v}_{a} & \delta\tilde{p}_{a} & \delta\tilde{c}_{a}^{\alpha} & \delta\tilde{c}_{a}^{\beta}\end{array}\right]\cdot\left[\begin{array}{cccc} \mathbf{K}_{ab}^{uu} & \mathbf{0} & \mathbf{0} & \mathbf{0}\\ \mathbf{k}_{ab}^{pu} & 0 & 0 & 0\\ \mathbf{k}_{ab}^{\alpha u} & 0 & 0 & 0\\ \mathbf{k}_{ab}^{\beta u} & 0 & 0 & 0 \end{array}\right]\cdot\left[\begin{array}{c} \Delta\mathbf{u}_{b}\\ \Delta\tilde{p}_{b}\\ \Delta\tilde{c}_{b}^{\alpha}\\ \Delta\tilde{c}_{b}^{\beta} \end{array}\right]\,,\label{eq346} \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}^{\alpha u} & =-\tilde{j}_{n}^{\alpha}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{eq347} \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})\).
Electric Potential and Partition Coefficient Derivatives
When the mixture is charged it is necessary to solve for the electric potential \(\psi\) using the electroneutrality condition in (2.9-4). This equation may be rewritted as a polynomial in \(\zeta\),
\[ \begin{equation} \sum\limits_{i=0}^{n}a_{i}\zeta^{i}\,,\label{eq348} \end{equation} \]
where
\[ \begin{equation} \zeta=\exp\left(-\frac{F_{c}\psi}{R\theta}\right)\,,\label{eq349} \end{equation} \]
and
\[ \begin{equation} a_{i}=\begin{cases} z^{\alpha}\hat{\kappa}^{\alpha}\tilde{c}^{\alpha} & i=z^{\alpha}-z^{\min}\\ c^{F} & i=-z^{\min} \end{cases}\,.\label{eq350} \end{equation} \]
Here, \(z^{\min}=\min_{\alpha}z^{\alpha}\) and the polynomial degress is \(n=z^{\max}-z^{\min}\) where \(z^{\max}=\max_{\alpha}z^{\alpha}\). Since more than one solute may carry the same charge \(z^{\alpha}\), the coefficients \(a_{i}\) should be evaluated from the summation of \(z^{\alpha}\hat{\kappa}^{\alpha}\tilde{c}^{\alpha}\) over all such solutes. Only real positive roots are valid, since \(\psi=-R\theta\left(\ln\zeta\right)/F_{c}\) according to \eqref{eq349}. Using Descartes' rule of signs, an inspection of the coefficients \(a_{i}\) shows tht there is only one sign change in the polynomial, regardless of the sign of \(c^{F}\), implying that there will always be only one positive root \(\zeta\), which must thus be real. Therefore, there cannot be any ambiguity in the calculation of \(\psi\), irrespective of the polynomial degree. Newton's method is used to solve for the positive real root when \(n>2\).
Using the above relations, it follows that \(\tilde{\kappa}^{\alpha}=\hat{\kappa}^{\alpha}\zeta^{z^{\alpha}}\). An examination of the equations resulting from the linearization of the internal virtual work shows that it is necessary to evaluate derivatives of \(\tilde{\kappa}^{\alpha}\) with respect to \(J\) and \(\tilde{c}^{\gamma}\), which are given by
\[ \begin{equation} \begin{aligned}\frac{\partial\tilde{\kappa}^{\alpha}}{\partial J} & =\frac{\partial\hat{\kappa}^{\alpha}}{\partial J}\zeta^{z^{\alpha}}+z^{\alpha}\tilde{\kappa}^{\alpha}\frac{1}{\zeta}\frac{\partial\zeta}{\partial J}\\ \frac{\partial\tilde{\kappa}^{\alpha}}{\partial\tilde{c}^{\gamma}} & =\frac{\partial\hat{\kappa}^{\alpha}}{\partial\tilde{c}^{\gamma}}\zeta^{z^{\alpha}}+z^{\alpha}\tilde{\kappa}^{\alpha}\frac{1}{\zeta}\frac{\partial\zeta}{\partial\tilde{c}^{\gamma}} \end{aligned} \,.\label{eq351} \end{equation} \]
In these expressions, the derivatives of \(\hat{\kappa}^{\alpha}\) are obtained from the user-defined constitutive relations for the solubility. Derivatives of \(\zeta\) may be evaluated by differentiating the electroneutrality condition to produce
\[ \begin{equation} \begin{aligned}\frac{1}{\zeta}\frac{\partial\zeta}{\partial J} & =-\frac{\frac{\partial c^{F}}{\partial J}+\sum\nolimits_{\beta}z^{\beta}\zeta^{z^{\beta}}\tilde{c}^{\beta}\frac{\partial\hat{\kappa}^{\beta}}{\partial J}}{\sum\nolimits_{\beta}\left(z^{\beta}\right)^{2}\tilde{\kappa}^{\beta}\tilde{c}^{\beta}}\\ \frac{1}{\zeta}\frac{\partial\zeta}{\partial\tilde{c}^{\gamma}} & =-\frac{z^{\gamma}\tilde{\kappa}^{\gamma}+\sum\nolimits_{\beta}z^{\beta}\zeta^{z^{\beta}}\tilde{c}^{\beta}\frac{\partial\hat{\kappa}^{\beta}}{\partial\tilde{c}^{\gamma}}}{\sum\nolimits_{\beta}\left(z^{\beta}\right)^{2}\tilde{\kappa}^{\beta}\tilde{c}^{\beta}} \end{aligned} \,.\label{eq352} \end{equation} \]
The derivative \(\partial c^{F}/\partial J\) may be evaluated from
\[ \begin{equation} c^{F}=\frac{1-\varphi_{r}^{s}}{J-\varphi_{r}^{s}}c_{r}^{F}\,,\label{eq353} \end{equation} \]
where \(\varphi_{r}^{s}\) is the referential solid volume fraction (volume of solid in current configuration per volume of the mixture in the reference configuration) and \(c_{r}^{F}\) is the referential fixed charge density (equivalent charge in current configuration per volume of the mixture in the reference configuration).
Chemical Reactions
Virtual Work and Linearization
The contribution to \(\delta W\) due to chemical reactions is given by \(\delta G\), where
\[ \divg\left(\mathbf{v}^{s}+\mathbf{w}\right)=\hat{\varphi}^{w}+\left(1-\varphi^{s}\right)\hat{\zeta}\sum_{\alpha}\nu^{\alpha}\mathcal{V}^{\alpha} \]
\[ \begin{equation} \delta G=\int_{b}\delta\tilde{p}\left[\hat{\varphi}^{w}+\left(1-\varphi^{s}\right)\hat{\zeta}\overline{\mathcal{V}}\right]\,dv+\sum\limits_{\iota}\nu^{\iota}\int_{b}\delta\tilde{c}^{\iota}\left(1-\varphi^{s}\right)\hat{\zeta}\,dv\,.\label{eq354} \end{equation} \]
The linearization of \(\delta G\) along a solid displacement increment \(\Delta\mathbf{u}\) is
\[ \begin{equation} \begin{aligned}D\delta G\left[\Delta\mathbf{u}\right] & =\int_{b}\delta\tilde{p}\,\left(\hat{\varphi}^{w}\,\divg\Delta\mathbf{u}+\hat{\boldsymbol{\varphi}}_{\varepsilon}^{w}:\Delta\boldsymbol{\varepsilon}\right)\,dv\\ & +\overline{\mathcal{V}}\int_{b}\delta\tilde{p}\left[\hat{\zeta}\,\divg\Delta\mathbf{u}+\left(J-\varphi_{r}^{s}\right)\hat{\boldsymbol{\zeta}}_{\varepsilon}:\Delta\boldsymbol{\varepsilon}\right]\,dv\\ & +\sum_{\iota}\nu^{\iota}\int_{b}\delta\tilde{c}^{\iota}\left[\hat{\zeta}\left(1-\frac{\partial\varphi_{r}^{s}}{\partial J}\right)\,\divg\Delta\mathbf{u}+\left(J-\varphi_{r}^{s}\right)\hat{\boldsymbol{\zeta}}_{\varepsilon}:\Delta\boldsymbol{\varepsilon}\right]\,dv \end{aligned} \,,\label{eq354b} \end{equation} \]
where
\[ \begin{equation} \hat{\boldsymbol{\varphi}}_{\varepsilon}^{w}=\mathbf{F}\cdot\frac{\partial\hat{\varphi}^{w}}{\partial\mathbf{E}}\cdot\mathbf{F}^{T}\,,\quad\hat{\boldsymbol{\zeta}}_{\varepsilon}=J^{-1}\mathbf{F}\cdot\frac{\partial\hat{\zeta}}{\partial\mathbf{E}}\cdot\mathbf{F}^{T}\,,\label{eq354c} \end{equation} \]
and
\[ \frac{\partial\varphi_{r}^{s}}{\partial J}=\frac{\sum_{\sigma}\frac{M^{\sigma}}{\rho_{T}^{\sigma}}c^{\sigma}\left(1+\sum_{\sigma}\frac{M^{\sigma}}{\rho_{T}^{\sigma}}c^{\sigma}\right)+\left(J-\varphi_{0}^{s}\right)\sum_{\sigma}\frac{M^{\sigma}}{\rho_{T}^{\sigma}}\frac{1}{\tilde{\kappa}^{\sigma}}\frac{\partial\tilde{\kappa}^{\sigma}}{\partial J}c^{\sigma}}{\left(1+\sum\limits_{\sigma}\frac{M^{\sigma}c^{\sigma}}{\rho_{T}^{\sigma}}\right)^{2}} \]
for solutes-as-SBMs, and
\[ \frac{\partial\varphi_{r}^{s}}{\partial J}=\frac{\sum_{\sigma}\frac{\rho_{r}^{\sigma}}{\rho_{T}^{\sigma}}\left(J\varphi^{w}+\sum_{\sigma}\frac{\rho_{r}^{\sigma}}{\rho_{T}^{\sigma}}\right)+J\varphi^{w}\left(J-\varphi_{0}^{s}\right)\sum_{\sigma}\frac{\rho_{r}^{\sigma}}{\rho_{T}^{\sigma}}\frac{1}{\tilde{\kappa}^{\sigma}}\frac{\partial\tilde{\kappa}^{\sigma}}{\partial J}}{\left(J\varphi^{w}+\sum\limits_{\sigma}\frac{\rho_{r}^{\sigma}}{\rho_{T}^{\sigma}}\right)^{2}} \]
for SBMs, as obtained from eq.(2.12-16). Currently, \(\hat{\zeta}\) is assumed to be independent of \(\tilde{p}\) in FEBio; it follows that the linearization along the effective fluid pressure increment \(\Delta\tilde{p}\) is
\[ \begin{equation} D\delta G\left[\Delta\tilde{p}\right]=\int_{b}\delta\tilde{p}\,\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{p}}\Delta p\,dv\,.\label{eq354d} \end{equation} \]
Finally, the linearization along a concentration increment \(\Delta\tilde{c}^{\iota}\) is
\[ \begin{equation} \begin{aligned}D\delta G\left[\Delta\tilde{c}^{\iota}\right] & =\int_{b}\delta\tilde{p}\,\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{c}^{\iota}}\Delta\tilde{c}^{\iota}\,dv\\ & +\overline{\mathcal{V}}\int_{b}\delta\tilde{p}\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}\Delta\tilde{c}^{\iota}\,dv\\ & +\sum_{\gamma}\nu^{\gamma}\int_{b}\delta\tilde{c}^{\gamma}\left(\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}-\frac{1}{J}\frac{\partial\varphi_{r}^{s}}{\partial\tilde{c}^{\iota}}\hat{\zeta}\right)\Delta\tilde{c}^{\iota}\,dv \end{aligned} \,.\label{eq354e} \end{equation} \]
where
\[ \frac{\partial\varphi_{r}^{s}}{\partial\tilde{c}^{\iota}}=\begin{cases} \frac{\left(J-\varphi_{r}^{s}\right)^{2}}{J-\varphi_{0}^{s}}\sum_{\sigma}\frac{M^{\sigma}}{\rho_{T}^{\sigma}}\left(\tilde{\kappa}^{\sigma}+\tilde{c}^{\sigma}\frac{\partial\tilde{\kappa}^{\sigma}}{\partial\tilde{c}^{\sigma}}\right) & \iota=\sigma\\ \frac{\left(J-\varphi_{r}^{s}\right)^{2}}{J-\varphi_{0}^{s}}\sum_{\sigma}\frac{M^{\alpha}}{\rho_{T}^{\sigma}}\tilde{c}^{\sigma}\frac{\partial\tilde{\kappa}^{\sigma}}{\partial\tilde{c}^{\iota}} & \iota\ne\sigma \end{cases} \]
for solutes-as-SBMs, as obtained from eq.(2.12-16).
The discretized form of these expressions is given by
\[ \begin{equation} \delta G=\sum_{a}\delta\tilde{p}_{a}r_{a}^{p}+\sum_{\gamma}\sum_{a}\delta\tilde{c}_{a}^{\gamma}r_{a}^{\gamma}\,,\label{eq354f} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}r_{a}^{p} & =\int_{b}N_{a}\,\left[\hat{\varphi}^{w}+\left(1-\varphi^{s}\right)\hat{\zeta}\bar{\mathcal{V}}\right]\,dv\\ r_{a}^{\gamma} & =\nu^{\gamma}\int_{b}N_{a}\left(1-\varphi^{s}\right)\hat{\zeta}\,dv \end{aligned} \,.\label{eq354g} \end{equation} \]
Similarly,
\[ \begin{equation} \begin{aligned}D\delta G\left[\Delta\mathbf{u}\right] & =\sum_{a}\delta\tilde{p}_{a}\sum_{b}\mathbf{k}_{ab}^{pu}\cdot\Delta\mathbf{u}_{b}\\ & +\sum_{\gamma}\sum_{a}\delta\tilde{c}_{a}^{\gamma}\sum_{b}\mathbf{k}_{ab}^{\gamma u}\cdot\Delta\mathbf{u}_{b} \end{aligned} \,,\label{eq354h} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}\mathbf{k}_{ab}^{pu} & =\int_{b}N_{a}\left(\hat{\varphi}^{w}\mathbf{I}+\hat{\boldsymbol{\varphi}}_{\varepsilon}^{w}\right)\cdot\grad N_{b}\,dv\\ & +\bar{\mathcal{V}}\int_{a}N_{a}\left[\hat{\zeta}\mathbf{I}+J\varphi^{w}\hat{\boldsymbol{\zeta}}_{\varepsilon}\right]\cdot\grad N_{b}\,dv\\ \mathbf{k}_{ab}^{\gamma u} & =\nu^{\gamma}\int_{a}N_{a}\left[\left(1-\frac{\partial\varphi_{r}^{s}}{\partial J}\right)\hat{\zeta}\mathbf{I}+J\varphi^{w}\hat{\boldsymbol{\zeta}}_{\varepsilon}\right]\cdot\grad N_{b}\,dv \end{aligned} \,.\label{eq354i} \end{equation} \]
Then,
\[ \begin{equation} D\delta G\left[\Delta\tilde{p}\right]=\sum_{a}\delta\tilde{p}_{a}\sum_{b}k_{ab}^{pp}\Delta\tilde{p}_{b}\,,\label{eq354j} \end{equation} \]
where
\[ \begin{equation} k_{ab}^{pp}=\int_{b}N_{a}\,\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{p}}N_{b}\,dv\,.\label{eq354k} \end{equation} \]
Finally,
\[ \begin{equation} \begin{aligned}D\delta G\left[\Delta\tilde{c}^{\iota}\right] & =\sum_{a}\delta\tilde{p}_{a}\sum_{b}k_{ab}^{p\iota}\Delta\tilde{c}_{b}^{\iota}\\ & +\sum_{\gamma}\sum_{a}\delta\tilde{c}_{a}^{\gamma}\sum_{b}k_{ab}^{\gamma\iota}\Delta\tilde{c}_{b}^{\iota} \end{aligned} \,,\label{eq354l} \end{equation} \]
where
\[ \begin{equation} \begin{aligned}k_{ab}^{p\iota} & =\int_{b}N_{a}N_{b}\frac{\partial\hat{\varphi}^{w}}{\partial\tilde{c}^{\iota}}\,dv+\bar{\mathcal{V}}\int_{b}N_{a}N_{b}\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}\,dv\\ k_{ab}^{\gamma\iota} & =\nu^{\gamma}\int_{b}N_{a}N_{b}\left(\left(1-\varphi^{s}\right)\frac{\partial\hat{\zeta}}{\partial\tilde{c}^{\iota}}-\frac{1}{J}\frac{\partial\varphi_{r}^{s}}{\partial\tilde{c}^{\iota}}\hat{\zeta}\right)\,dv \end{aligned} \,.\label{eq354m} \end{equation} \]
Updating Solid-Bound Molecule Concentrations
The solid-bound molecule concentrations \(\rho_{r}^{\sigma}\) are evaluated at integration points of each element; they do not represent nodal degrees of freedom. The values of \(\rho_{r}^{\sigma}\) are updated at the end of each iteration in the solution of the nonlinear equations for the nodal degrees of freedom, using trapezoidal integration on \(\hat{\rho}_{r}^{\sigma}\) in (2.12-6). According to (2.12-24) and (2.12-26), we have \(\hat{\rho}_{r}^{\sigma}=\left(J-\varphi_{r}^{s}\right)M^{\sigma}\nu^{\sigma}\hat{\zeta}\), which is evaluated as the average of values at \(t_{n}\) and \(t_{n+1}\), then
\[ \begin{equation} \left(\rho_{r}^{\sigma}\right)_{n+1}=\left(\rho_{r}^{\sigma}\right)_{n}+\left(\hat{\rho}_{r}^{\sigma}\right)_{n+\frac{1}{2}}\Delta t\label{eq:sbm-incremental-eq} \end{equation} \]
where \(\Delta t=t_{n+1}-t_{n}\).