To linearize the finite element equations, the directional derivative of the virtual work in equation (3.1-1) must be calculated. In an iterative procedure, the quantity $\boldsymbol{\phi}$ will be approximated by a trial solution $\boldsymbol{\phi}_{k}$ . Linearization of the virtual work equation around this trial solution gives $$\begin{equation} \delta W\left(\boldsymbol{\phi}_{k},\delta\mathbf{v}\right)+D\delta W\left(\boldsymbol{\phi}_{k},\delta\mathbf{v}\right)\left[\mathbf{u}\right]=0\,.\label{eq173} \end{equation}$$ The directional derivative of the virtual work will eventually lead to the definition of the stiffness matrix. In order to proceed, it is convenient to split the virtual work into an internal and external virtual work component: $$\begin{equation} D\delta W\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right]=D\delta W_{int}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right]-D\delta W_{ext}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right],\label{eq174} \end{equation}$$ where $$\begin{equation} \delta W_{int}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)=\int\limits _{v}\boldsymbol{\sigma}:\delta\mathbf{d}\,dv\,,\label{eq175} \end{equation}$$ and $$\begin{equation} \delta W_{ext}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)=\int\limits _{v}\mathbf{f}\cdot\delta\mathbf{v}\,dv+\int\limits _{\partial v}\mathbf{t}\cdot\delta\mathbf{v}\,da\,.\label{eq176} \end{equation}$$ The result is listed here without details of the derivation – see [23] for details. The linearization of the internal virtual work is given by $$\begin{equation} D\delta W_{int}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right]=\int\limits _{v}\delta\mathbf{d}:\boldsymbol{\mathcal{C}}:\boldsymbol{\varepsilon}\,dv+\int\limits _{v}\boldsymbol{\sigma}:\left[\left(\nabla\mathbf{u}\right)^{T}\cdot\nabla\delta\mathbf{v}\right]\,dv\,.\label{eq177} \end{equation}$$ Notice that this equation is symmetric in $\delta\mathbf{v}$ and $\mathbf{u}$ . This symmetry will, upon discretization, yield a symmetric tangent matrix.

The external virtual work has contributions from both body forces and surface tractions. The precise form of the linearized external virtual work depends on the form of these forces. For surface tractions, normal pressure forces may be represented in FEBio. The linearized external work for this type of traction is given by $$\begin{equation} \begin{aligned}D\delta W_{ext}^{p}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right] & =\frac{1}{2}\int\limits _{A_{\xi}}p\frac{\partial\mathbf{x}}{\partial\xi}\cdot\left[\left(\frac{\partial\mathbf{u}}{\partial\eta}\times\delta\mathbf{v}\right)+\left(\frac{\partial\delta\mathbf{v}}{\partial\eta}\times\mathbf{u}\right)\right]d\xi\,d\eta\\ & -\frac{1}{2}\int\limits _{A_{\xi}}p\frac{\partial\mathbf{x}}{\partial\eta}\cdot\left[\left(\frac{\partial\mathbf{u}}{\partial\xi}\times\delta\mathbf{v}\right)+\left(\frac{\partial\delta\mathbf{v}}{\partial\xi}\times\mathbf{u}\right)\right]d\xi\,d\eta\,. \end{aligned} \label{eq178} \end{equation}$$ Discretization of this equation will also lead to a symmetric component of the tangent matrix.

FEBio currently supports gravity as a body force, $\mathbf{f}=\rho\mathbf{g}$ . Since this force is independent of the geometry, the contribution to the linearized external work is zero. Another type of body force implemented in FEBio is the centrifugal force. For a body rotating with a constant angular speed $\omega$ , about an axis passing through the point $\mathbf{c}$ and directed along the unit vector $\mathbf{n}$ , the body force is given by $\mathbf{f}=\rho\omega^{2}\mathbf{r}$ , where $\mathbf{r}$ is the vector distance from a point $\mathbf{x}$ to the axis of rotation, $$\begin{equation} \mathbf{r}=\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)\cdot\left(\mathbf{x}-\mathbf{c}\right)\,.\label{eq179} \end{equation}$$

The resulting linearized external work is given by $$\begin{equation} D\delta W_{ext}^{f}\left(\boldsymbol{\phi},\delta\mathbf{v}\right)\left[\mathbf{u}\right]=\int\limits _{v}\rho\omega^{2}\delta\mathbf{v}\cdot\left(\mathbf{I}-\mathbf{n}\otimes\mathbf{n}\right)\cdot\mathbf{u}\,dv,\label{eq180} \end{equation}$$ which produces a symmetric expression that will yield a symmetric matrix.