$\newcommand{\Dev}{\operatorname{Dev}}$ $\newcommand{\dev}{\operatorname{dev}}$ $\newcommand{\grad}{\operatorname{grad}}$ $\newcommand{\Grad}{\operatorname{Grad}}$ $\newcommand{\divg}{\operatorname{div}}$ $\newcommand{\Ei}{\operatorname{Ei}}$ $\newcommand{\tr}{\operatorname{tr}}$ $\newcommand{\Divg}{\operatorname{Div}}$ $\newcommand{\cay}{\operatorname{cay}}$ $\newcommand{\rot}{\operatorname{rot}}$

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.7.1-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 3.3.1↑-3.3.2↑ 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.7.1-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.