6.2 Elastodynamics
Governing Equations
The linear momentum balance for elastodynamics is
\[ \begin{equation} \rho\mathbf{a}=\divg\boldsymbol{\sigma}+\rho\mathbf{b}\,,\label{eq:edy-solid-momentum} \end{equation} \]
where \(\rho\) is the density, \(\mathbf{a}\) is the acceleration, \(\boldsymbol{\sigma}\) is the Cauchy stress, and \(\mathbf{b}\) is the body force per mass. The angular momentum balance is satisfied by letting \(\boldsymbol{\sigma}^{T}=\boldsymbol{\sigma}\). The integrated form of the mass balance equations yields
\[ \begin{equation} \rho=\frac{\rho_{r}}{J}\,,\label{eq:edy-mass-balance} \end{equation} \]
where \(\rho_{r}\) is the density in the reference configuration and \(J=\det\mathbf{F}\), where \(\mathbf{F}\) is the deformation gradient. The acceleration is given by the material time derivative of the velocity \(\mathbf{v}\), evaluated either in a spatial or a material frame,
\[ \begin{equation} \mathbf{a}=\dot{\mathbf{v}}\label{eq:edy-acceleration} \end{equation} \]
Virtual Work
The virtual work for the domain \(b\) is given by
\[ \begin{equation} \delta W=\int_{b}\delta\mathbf{v}\cdot\left(\divg\boldsymbol{\sigma}+\rho\left(\mathbf{b}-\mathbf{a}\right)\right)\,dv\,,\label{eq:edy-virtual-work} \end{equation} \]
where \(\delta\mathbf{v}\) is the virtual velocity. Using the divergence theorem, this virtual work may be expressed as the difference \(\delta W=\delta W_{ext}-\delta W_{int}\) between external virtual work \(\delta W_{ext}\) and internal virtual work \(\delta W_{int}\), where
\[ \begin{equation} \begin{aligned}\delta W_{int} & =\int_{b}\boldsymbol{\sigma}:\grad\delta\mathbf{v}\,dv+\int_{b}\delta\mathbf{v}\cdot\rho\mathbf{a}\,dv\,,\\ \delta W_{ext} & =\int_{\partial b}\delta\mathbf{v}\cdot\mathbf{t}\,da+\int_{b}\delta\mathbf{v}\cdot\rho\mathbf{b}\,dv\,, \end{aligned} \label{eq:edy-int-ext-virtual-work} \end{equation} \]
where \(\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n}\) is the traction on the boundary \(\partial b\).
Generalized \(\alpha-\)Method for Elastodynamics
In the generalized \(\alpha-\)method, we evaluate displacements and velocities at the intermediate time \(t_{n+\alpha_{f}}=t_{n}+\alpha_{f}\left(t_{n+1}-t_{n}\right)\), where \(\alpha_{f}\) is a user-defined parameter (\(0<\alpha_{f}\le1\)), such that
\[ \begin{equation} \begin{aligned}\boldsymbol{\chi}_{n+\alpha_{f}} & =\left(1-\alpha_{f}\right)\boldsymbol{\chi}_{n}+\alpha_{f}\boldsymbol{\chi}_{n+1}\,,\\ \mathbf{u}_{n+\alpha_{f}} & =\left(1-\alpha_{f}\right)\mathbf{u}_{n}+\alpha_{f}\mathbf{u}_{n+1}\,,\\ \mathbf{v}_{n+\alpha_{f}} & =\left(1-\alpha_{f}\right)\mathbf{v}_{n}+\alpha_{f}\mathbf{v}_{n+1}\,, \end{aligned} \label{eq:edy-motion-interpolation} \end{equation} \]
where \(\boldsymbol{\chi}\) is the motion and \(\mathbf{u}\) is the displacement. In particular, it follows that the deformation gradient and its determinant are given at the intermediate time by
\[ \begin{equation} \mathbf{F}_{n+\alpha_{f}}=\frac{\partial\boldsymbol{\chi}_{n+\alpha_{f}}}{\partial\mathbf{X}}=\left(1-\alpha_{f}\right)\mathbf{F}_{n}+\alpha_{f}\mathbf{F}_{n+1}\,,\label{eq:edy-mp-defgrad} \end{equation} \]
and
\[ \begin{equation} J_{n+\alpha_{f}}=\det\mathbf{F}_{n+\alpha_{f}}\,.\label{eq:edy-mp-detJ} \end{equation} \]
The material time derivative of \(J_{n+\alpha_{f}}\), and the velocity gradient \(\mathbf{L}_{n+\alpha_{f}}\) are normally evaluated as
\[ \begin{equation} \dot{J}_{n+\alpha_{f}}=J_{n+\alpha_{f}}\mathbf{F}_{n+\alpha_{f}}^{-T}:\Grad\mathbf{v}_{n+\alpha_{f}}\,,\label{eq:edy-Jdot-exact} \end{equation} \]
and
\[ \begin{equation} \mathbf{L}_{n+\alpha_{f}}=\Grad\mathbf{v}_{n+\alpha_{f}}\cdot\mathbf{F}_{n+\alpha_{f}}^{-1}\,.\label{eq:edy-L-exact} \end{equation} \]
In practice however, we get better numerical results when using
\[ \begin{equation} \dot{J}_{n+\alpha_{f}}=\frac{J_{n+1}-J_{n}}{\Delta t}\,,\label{eq:edy-Jdot-approx} \end{equation} \]
and
\[ \begin{equation} \mathbf{L}_{n+\alpha_{f}}=\frac{\mathbf{F}_{n+1}-\mathbf{F}_{n}}{\Delta t}\cdot\mathbf{F}_{n+\alpha_{f}}^{-1}\,.\label{eq:edy-L-approx} \end{equation} \]
According to the generalized\(-\alpha\) method, we evaluate the velocity derivative at a different intermediate time \(t_{n+\alpha_{m}}=t_{n}+\alpha_{m}\left(t_{n+1}-t_{n}\right)\), such that
\[ \begin{equation} \dot{\mathbf{v}}_{n+\alpha_{m}}=\left(1-\alpha_{m}\right)\dot{\mathbf{v}}_{n}+\alpha_{m}\dot{\mathbf{v}}_{n+1}\,.\label{eq:edy-vdot-interpolation} \end{equation} \]
Since elastodynamics represent a second-order system of equations in time, the parameters \(\alpha_{f}\) and \(\alpha_{m}\) are evaluated from a single parameter \(\rho_{\infty}\) using
\[ \begin{equation} \alpha_{f}=\frac{1}{1+\rho_{\infty}}\,,\quad\alpha_{m}=\frac{2-\rho_{\infty}}{1+\rho_{\infty}}\,,\label{eq:ga-alphas-2-1} \end{equation} \]
where \(0\le\rho_{\infty}\le1\). This parameter is the spectral radius for an infinite time step, which controls the amount of damping of high frequencies; a value of zero produces the greatest amount of damping, anihilating the highest frequency in one step, whereas a value of one preserves the highest frequency.
To complete the integration scheme , we evaluate
\[ \begin{equation} \begin{aligned}\beta & =\frac{1}{4}\left(1+\alpha_{m}-\alpha_{f}\right)^{2}\\ \gamma & =\frac{1}{2}+\alpha_{m}-\alpha_{f} \end{aligned} \,,\label{eq:edy-Newmark} \end{equation} \]
then we use the Newmark integration formulas (Section Newmark Integration),
\[ \begin{equation} \begin{aligned}\mathbf{v}_{n+1} & =\mathbf{v}_{n}+\Delta t\left[\left(1-\gamma\right)\dot{\mathbf{v}}_{n}+\gamma\dot{\mathbf{v}}_{n+1}\right]\\ \mathbf{u}_{n+1} & =\mathbf{u}_{n}+\Delta t\mathbf{v}_{n}+\frac{\Delta t^{2}}{2}\left[\left(1-2\beta\right)\dot{\mathbf{v}}_{n}+2\beta\dot{\mathbf{v}}_{n+1}\right]\\ \dot{\mathbf{v}}_{n+1} & =\frac{1}{\beta\Delta t}\left(\frac{\mathbf{u}_{n+1}-\mathbf{u}_{n}}{\Delta t}-\mathbf{v}_{n}\right)+\left(1-\frac{1}{2\beta}\right)\dot{\mathbf{v}}_{n} \end{aligned} \,.\label{eq:edy-current-time} \end{equation} \]
At the start of each time step, we initialize the variables as follows:
\[ \begin{equation} \begin{aligned}\mathbf{u}_{n+1} & =\mathbf{u}_{n}\\ \dot{\mathbf{v}}_{n+1} & =\left(1-\frac{1}{2\beta}\right)\dot{\mathbf{v}}_{n}-\frac{1}{\beta\Delta t}\mathbf{v}_{n}\\ \mathbf{v}_{n+1} & =\left(1-\frac{\gamma}{\beta}\right)\mathbf{v}_{n}+\Delta t\left(1-\frac{\gamma}{2\beta}\right)\dot{\mathbf{v}}_{n} \end{aligned} \,.\label{eq:edy-initializations} \end{equation} \]
Linearization
The solution of the nonlinear equation \(\delta W=0\) is obtained by linearizing this relation as
\[ \begin{equation} \delta W+D\delta W\left[\Delta\mathbf{u}\right]\approx0\,,\label{eq:edy-linearized-virtual-work} \end{equation} \]
where the operator \(D\delta W\left[\cdot\right]\) represents the directional derivative of \(\delta W\) at \(\mathbf{u}\) along an increment \(\Delta\mathbf{u}\) of \(\mathbf{u}_{n+1}\) . According to the generalized\(-\alpha\) method , the virtual work is evaluated using intermediate time step values, at \(t_{n+\alpha_{f}}\) for all parameters except \(\dot{\mathbf{v}}\), which is evaluated at \(t_{n+\alpha_{m}}\). It follows from these definitions that the linearizations of critical variables are given by
\[ \begin{equation} \begin{aligned}D\mathbf{u}\left[\Delta\mathbf{u}\right] & =\alpha_{f}\Delta\mathbf{u}\\ D\mathbf{F}\left[\Delta\mathbf{u}\right] & =\alpha_{f}\Grad\Delta\mathbf{u}\\ DJ\left[\Delta\mathbf{u}\right] & =\alpha_{f}J\left(\divg\Delta\mathbf{u}\right)\\ D\dot{J}\left[\Delta\mathbf{u}\right] & =D\left(J\mathbf{F}^{-T}:\Grad\mathbf{v}\right)\left[\Delta\mathbf{u}\right]\\ & =\alpha_{f}J\left[\left(\divg\mathbf{v}+\frac{\gamma}{\beta\Delta t}\right)\left(\divg\Delta\mathbf{u}\right)-\left(\grad\Delta\mathbf{u}\right)^{T}:\mathbf{L}\right]\\ D\mathbf{v}\left[\Delta\mathbf{u}\right] & =\frac{\alpha_{f}\gamma}{\beta\Delta t}\Delta\mathbf{u}\\ D\dot{\mathbf{v}}\left[\Delta\mathbf{u}\right] & =\frac{\alpha_{m}}{\beta\Delta t^{2}}\Delta\mathbf{u} \end{aligned} \,.\label{eq:edy-linearizations} \end{equation} \]
To linearize the virtual work, we need to express the integrals appearing in \(\delta W_{int}\) and \(\delta W_{ext}\) over the material frame of the finite element solid domain.
Internal Work
The first term in the internal work becomes
\[ \begin{equation} \begin{aligned}\int_{b}\boldsymbol{\sigma}:\grad\delta\mathbf{v}\,dv & =\int_{B}\mathbf{F}\cdot\mathbf{S}:\Grad\delta\mathbf{v}\,dV\end{aligned} \,,\label{eq:work-material-1-1} \end{equation} \]
where \(\mathbf{S}=J\cdot\mathbf{F}^{-1}\cdot\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\) is the second Piola-Kirchhoff stress for the solid material. In general, \(\boldsymbol{\sigma}\) (and thus, \(\mathbf{S}\)) is only a function of the solid strain, such as the right Cauchy-Green tensor \(\mathbf{C}=\mathbf{F}^{T}\cdot\mathbf{F}\) or the Green-Lagrange strain \(\mathbf{E}=\left(\mathbf{C}-\mathbf{I}\right)/2\).
\[ \begin{equation} D\mathbf{E}\left[\Delta\mathbf{u}\right]=\frac{\alpha_{f}}{2}\left(\Grad^{T}\Delta\mathbf{u}\cdot\mathbf{F}+\mathbf{F}^{T}\cdot\Grad\Delta\mathbf{u}\right)\,.\label{eq:edy-E-linearization} \end{equation} \]
Therefore, following the standard approach in solid mechanics, the linearization of \(\mathbf{S}\) is
\[ \begin{equation} \begin{aligned}D\mathbf{S}\left[\Delta\mathbf{u}\right] & =\frac{\partial\mathbf{S}}{\partial\mathbf{E}}:D\mathbf{E}\left[\Delta\mathbf{u}\right]\\ & =\alpha_{f}\boldsymbol{\mathbb{C}}:\frac{1}{2}\left(\Grad^{T}\Delta\mathbf{u}\cdot\mathbf{F}+\mathbf{F}^{T}\cdot\Grad\Delta\mathbf{u}\right)\\ & =\alpha_{f}\boldsymbol{\mathbb{C}}:\left(\mathbf{F}^{T}\oslash\mathbf{F}^{T}\right):\Delta\boldsymbol{\varepsilon} \end{aligned} \,,\label{eq:edy-solid-stress-linearization} \end{equation} \]
where \(\boldsymbol{\mathbb{C}}\) is the material elasticity tensor. Now, the linearization of the first term in \(\delta W_{int}\) is
\[ \begin{equation} \boxed{\begin{aligned} & D\left(\int_{B}\mathbf{F}\cdot\mathbf{S}:\Grad\delta\mathbf{v}\,dV\right)\left[\Delta\mathbf{u}\right]\\ & =\int_{v}\alpha_{f}\left(\grad\delta\mathbf{v}:\grad\Delta\mathbf{u}\cdot\boldsymbol{\sigma}+\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}:\grad\Delta\mathbf{u}\right)\,dv \end{aligned} }\,.\label{eq:edy-Dint-1} \end{equation} \]
where \(\boldsymbol{\mathcal{C}}\) is the spatial elasticity tensor. Similarly, the second term in \(\delta W_{int}\) produces
\[ \begin{equation} \boxed{D\left(\int_{B}\delta\mathbf{v}\cdot\rho_{r}\mathbf{a}\,dV\right)\left[\Delta\mathbf{u}\right]=\int_{b}\delta\mathbf{v}\cdot\frac{\alpha_{m}}{\beta\Delta t^{2}}\rho\Delta\mathbf{u}\,dv}\,.\label{eq:edy-Dint-2} \end{equation} \]
External Work
The linearization of the body force term is
\[ \begin{equation} D\left(\int_{B}\delta\mathbf{v}\cdot\rho_{r}\mathbf{b}\,dV\right)\left[\Delta\mathbf{u}\right]=\int_{b}\delta\mathbf{v}\cdot\alpha_{f}\rho\grad\mathbf{b}\cdot\Delta\mathbf{u}\,dv\,.\label{eq:edy-Dext} \end{equation} \]
The linearization of the traction force term is
\[ \begin{equation} \begin{aligned} & D\left(\int_{\Gamma_{\eta}}\delta\mathbf{v}\cdot\mathbf{t}\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\,d\eta^{1}d\eta^{2}\right)\left[\Delta\mathbf{u}\right]\\ & =\int_{\Gamma_{\eta}}\alpha_{f}\delta\mathbf{v}\cdot\left(\mathbf{t}\otimes\mathbf{n}\right)\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)\,d\eta^{1}d\eta^{2}\\ & +\int_{\Gamma_{\eta}}\delta\mathbf{v}\cdot D\mathbf{t}\left[\Delta\mathbf{u}\right]\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\,d\eta^{1}d\eta^{2} \end{aligned} \,.\label{eq:edy-Dext-final} \end{equation} \]
Note that \(D\mathbf{t}\left[\Delta\mathbf{u}\right]\) depends on the nature of the surface traction. For a prescribed traction we have \(D\mathbf{t}\left[\Delta\mathbf{u}\right]=\mathbf{0}\). A contact analysis needs more elaborate derivations (not yet implemented as of FEBio 2.7). In the above expression we used
\[ \begin{equation} \begin{aligned}D\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|\left[\Delta\mathbf{u}\right] & =\mathbf{n}\cdot D\left(\mathbf{g}_{1}\times\mathbf{g}_{2}\right)\left[\Delta\mathbf{u}\right]\\ & =\alpha_{f}\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) \end{aligned} \label{eq:edy-Dext-DA} \end{equation} \]
where the unit outward normal is evaluated as
\[ \begin{equation} \mathbf{n}=\frac{\mathbf{g}_{1}\times\mathbf{g}_{2}}{\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,.\label{eq:edy-unit-normal} \end{equation} \]
Discretization
We use the following interpolations:
\[ \begin{equation} \begin{aligned}\delta\mathbf{v} & =\sum_{a}N_{a}\delta\mathbf{v}_{a} & \Delta\mathbf{u} & =\sum_{b}N_{b}\Delta\mathbf{u}_{b}\\ \grad\delta\mathbf{v} & =\sum_{a}\delta\mathbf{v}_{a}\otimes\grad N_{a} & \grad\Delta\mathbf{u} & =\sum_{b}\Delta\mathbf{u}_{b}\otimes\grad N_{b}\\ \divg\delta\mathbf{v} & =\sum_{a}\delta\mathbf{v}_{a}\cdot\grad N_{a} & \divg\Delta\mathbf{u} & =\sum_{b}\Delta\mathbf{u}_{b}\cdot\grad N_{b} \end{aligned} \,,\label{eq:interpolations-1} \end{equation} \]
where \(N_{a}\left(\eta^{1},\eta^{2},\eta^{3}\right)\) are shape functions of the element parametric coordinates \(\left(\eta^{1},\eta^{2},\eta^{3}\right)\). Note that the \(\grad\equiv\frac{\partial}{\partial\mathbf{x}}\) operator should be evaluated at \(t_{n+\alpha_{f}}\), using \(\mathbf{x}_{n+\alpha_{f}}\). For example, in the case of a scalar function \(f\),
\[ \begin{aligned}\grad f & =\frac{\partial f}{\partial\mathbf{x}_{n+\alpha_{f}}}=\frac{\partial f}{\partial\eta^{i}}\mathbf{g}_{n+\alpha_{f}}^{i}\\ \mathbf{g}_{n+\alpha_{f}}^{i} & =\frac{\partial\eta^{i}}{\partial\mathbf{x}_{n+\alpha_{f}}} \end{aligned} \,, \]
where the contravariant basis vectors \(\mathbf{g}_{n+\alpha_{f}}^{i}\) may be evaluated from the covariant basis vectors
\[ \mathbf{g}_{i}^{n+\alpha_{f}}=\frac{\partial\mathbf{x}_{n+\alpha_{f}}}{\partial\eta^{i}}=\left(1-\alpha_{f}\right)\frac{\partial\mathbf{x}_{n}}{\partial\eta^{i}}+\alpha_{f}\frac{\partial\mathbf{x}_{n+1}}{\partial\eta^{i}} \]
using \(\mathbf{g}_{i}^{n+\alpha_{f}}\cdot\mathbf{g}_{n+\alpha_{f}}^{j}=\delta_{i}^{j}\).
The discretization of the internal work produces
\[ \begin{equation} \delta W_{int}=\sum_{a}\delta\mathbf{v}_{a}\cdot\int_{b}\left(\mathbf{f}_{a}^{u}+\mathbf{f}_{a}^{\rho}\right)\,dv\,,\label{eq:edy-discretized-internal-work} \end{equation} \]
where
\[ \begin{equation} \boxed{\begin{aligned}\mathbf{f}_{a}^{u} & =\boldsymbol{\sigma}\cdot\grad N_{a}\\ \mathbf{f}_{a}^{\rho} & =N_{a}\rho\mathbf{a} \end{aligned} }\,.\label{eq:edy-discretized-residuals} \end{equation} \]
The discretization of the stress and elasticity terms in the internal work is
\[ \begin{aligned} & \int_{v}\alpha_{f}\left(\grad\delta\mathbf{v}:\grad\Delta\mathbf{u}\cdot\boldsymbol{\sigma}+\grad\delta\mathbf{v}:\boldsymbol{\mathcal{C}}:\grad\Delta\mathbf{u}\right)\,dv\\ & =\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\int_{v}\mathbf{K}_{ab}\,dv\cdot\Delta\mathbf{u}_{b} \end{aligned} \,, \]
where
\[ \begin{equation} \boxed{\mathbf{K}_{ab}=\alpha_{f}\left(\left(\grad N_{a}\cdot\boldsymbol{\sigma}\cdot\grad N_{b}\right)\mathbf{I}+\grad N_{a}\cdot\boldsymbol{\mathcal{C}}\cdot\grad N_{b}\right)}\,.\label{eq:edy-Wint-Kab} \end{equation} \]
The discretization of the mass term in the internal work is
\[ \int_{b}\delta\mathbf{v}\cdot\frac{\alpha_{m}}{\beta\Delta t^{2}}\rho\Delta\mathbf{u}\,dv=\sum_{a}\delta\mathbf{v}_{a}\cdot\sum_{b}\int_{b}\mathbf{M}_{ab}\,dv\cdot\Delta\mathbf{u}_{b}\,, \]
where
\[ \begin{equation} \boxed{\mathbf{M}_{ab}=\frac{\alpha_{m}}{\beta\Delta t^{2}}\rho N_{a}N_{b}\mathbf{I}}\,.\label{eq:edy-Wint-Mab} \end{equation} \]
For the external work of body forces,
\[ \int_{b}\delta\mathbf{v}\cdot\rho\mathbf{b}\,dv=\sum_{a}\delta\mathbf{v}_{a}\cdot\int_{b}\mathbf{f}_{a}^{\mathbf{b}}\,dv \]
where
\[ \begin{equation} \boxed{\mathbf{f}_{a}^{\mathbf{b}}=N_{a}\rho^{s}\mathbf{b}}\,,\label{eq:edy-Wext-fb} \end{equation} \]
and
\[ \int_{b}\delta\mathbf{v}^{s}\cdot\alpha_{f}\rho^{s}\grad\mathbf{b}\cdot\Delta\mathbf{u}\,dv=\sum_{a}\delta\mathbf{v}_{a}^{s}\cdot\sum_{b}\int_{b}\mathbf{K}_{ab}^{\mathbf{b}}\,dv\cdot\Delta\mathbf{u}_{b} \]
where
\[ \begin{equation} \boxed{\mathbf{K}_{ab}^{\mathbf{b}}=\alpha_{f}N_{a}N_{b}\rho^{s}\grad\mathbf{b}}\,.\label{eq:edy-Wext-Kb} \end{equation} \]
For prescribed tractions,
\[ \int_{\partial b}\delta\mathbf{v}^{s}\cdot\mathbf{t}^{s}\,da=\sum_{a}\delta\mathbf{v}_{a}^{s}\cdot\int_{\Gamma_{\eta}}\mathbf{f}_{a}^{t}\,d\eta^{1}d\eta^{2} \]
where
\[ \begin{equation} \boxed{\mathbf{f}_{a}^{t}=N_{a}\mathbf{t}^{s}\,\left|\mathbf{g}_{1}\times\mathbf{g}_{2}\right|}\,,\label{eq:edy-Wext-ft} \end{equation} \]
and
\[ \begin{aligned} & \int_{\Gamma_{\eta}}\alpha_{f}\delta\mathbf{v}^{s}\cdot\left(\mathbf{t}^{s}\otimes\mathbf{n}\right)\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)\,d\eta^{1}d\eta^{2}\\ & =\sum_{a}\delta\mathbf{v}_{a}^{s}\cdot\sum_{b}\int_{\Gamma_{\eta}}\mathbf{K}_{ab}^{t}\,d\eta^{1}d\eta^{2}\cdot\Delta\mathbf{u}_{b} \end{aligned} \]
where
\[ \begin{equation} \boxed{\mathbf{K}_{ab}^{t}=\alpha_{f}N_{a}\left(\mathbf{t}^{s}\otimes\mathbf{n}\right)\cdot\left(\frac{\partial N_{b}}{\partial\eta^{2}}\hat{\mathbf{g}}_{1}-\frac{\partial N_{b}}{\partial\eta^{1}}\hat{\mathbf{g}}_{2}\right)}\,,\label{eq:edy-Wext-Kt} \end{equation} \]
where \(\hat{\mathbf{g}}\) is the skew-symmetric tensor whose dual vector is \(\mathbf{g}\).
Energy-Momentum Conservation Scheme
The time discretization scheme may be selected in a manner that enforces linear and angular momentum, and energy conservation over consecutive time steps \(t_{n}\) and \(t_{n+1}\), when boundary conditions and external loads are time-independent. Based on the prior literature , this momentum and energy conservation may be achieved by using the midpoint rule (\(\rho_{\infty}=1\), leading to \(\alpha_{f}=\alpha_{m}=\frac{1}{2}\)), and evaluating the virtual work at \(t_{n+\frac{1}{2}}\). However, since the virtual work strictly enforces momentum balance only, there is no guarantee that energy conservation will be satisfied as a result of time discretization. Therefore, we need to enforce a specific scheme to satisfy energy balance.
Energy Balance
For an elastic solid, in the absence of heat exchanges (i.e., in elastodynamics), the equation of energy balance reduces
\[ \begin{equation} \rho\dot{\varepsilon}=\boldsymbol{\sigma}:\mathbf{D}\,,\label{eq:edy-energy-balance} \end{equation} \]
where \(\varepsilon\) is the specific internal energy and \(\mathbf{D}\) is the rate of deformation tensor. Recall that \(\varepsilon=\psi+\theta\eta\), where \(\psi\) is the specific free energy, \(\theta\) is the absolute temperature and \(\eta\) is the specific entropy. Since \(\eta=0\) in elasticity (due to the temperature remaining constant), the above energy balance may be combined with the mass balance \eqref{eq:edy-mass-balance} as
\[ \begin{equation} \rho\dot{\psi}=\frac{\rho_{r}}{J}\dot{\psi}=\boldsymbol{\sigma}:\mathbf{D}\,,\label{eq:edy-energy-r1} \end{equation} \]
or
\[ \begin{equation} \dot{\Psi}_{r}=J\boldsymbol{\sigma}:\mathbf{D}\,,\label{eq:edy-energy-r2} \end{equation} \]
where \(\Psi_{r}=\rho_{r}\psi\) is the free energy density (per volume of the material in the reference configuration).
In our time integration scheme, to satisfy energy balance, this equation needs to be evaluated at \(t_{n+\alpha_{f}}\), thus
\[ \begin{equation} \left(\dot{\Psi}_{r}\right)_{n+\alpha_{f}}=J_{n+\alpha_{f}}\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\,.\label{eq:edy-mp-energy} \end{equation} \]
However, the solution for \(\boldsymbol{\sigma}_{n+\alpha_{f}}\equiv\boldsymbol{\sigma}\left(\mathbf{F}_{n+\alpha_{f}}\right)\) obtained from the momentum balance may not necessarity satisfy this equation. Thus, to satisfy energy balance over consecutive time steps, we want to evaluate an effective stress \(\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}\) such that
\[ \begin{equation} J_{n+\alpha_{f}}\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}=\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{\Delta t}\,.\label{eq:edy-energy-discretized} \end{equation} \]
To find a solution for \(\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}\), we follow the procedure of Gonzalez and let
\[ \begin{equation} \tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}+f\mathbf{D}_{n+\alpha_{f}}\,,\label{eq:edy-eff-stress-model} \end{equation} \]
where \(f\) is some scalar function to be determined. Substituting this relation, \eqref{eq:edy-eff-stress-model}, into the previous equation, \eqref{eq:edy-energy-discretized}, produces
\[ \begin{equation} f=\left(\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{J_{n+\alpha_{f}}^{s}\Delta t}-\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\right)\frac{1}{\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}}\,.\label{eq:edy-eff-stress-f} \end{equation} \]
Hence, the equation for an effective stress needed to satisfy energy balance between consecutive time steps is
\[ \begin{equation} \boxed{\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}+\left(\frac{\left(\Psi_{r}\right)_{n+1}-\left(\Psi_{r}\right)_{n}}{J_{n+\alpha_{f}}\Delta t}-\boldsymbol{\sigma}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}\right)\frac{\mathbf{D}_{n+\alpha_{f}}}{\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}}}\,.\label{eq:edy-eff-stress-final} \end{equation} \]
In the limit when \(\mathbf{D}_{n+\alpha_{f}}:\mathbf{D}_{n+\alpha_{f}}=0\), we use \(\tilde{\boldsymbol{\sigma}}_{n+\alpha_{f}}=\boldsymbol{\sigma}_{n+\alpha_{f}}\). Recall that this scheme produces conservation of linear and angular momentum and total energy only with \(\rho_{\infty}=1\), or equivalently, \(\alpha_{f}=\alpha_{m}=\frac{1}{2}\), \(\beta=\frac{1}{2}\) and \(\gamma=1\). Therefore, this effective stress calculation is only applied when the user employs \(\rho_{\infty}=1\).