  Computational biomechanics plays an important role in biomedical engineering:
using modeling to understand pathophysiology, treatment and device design.
While experimental evidence indicates that the mechanical response of most
tissues is viscoelasticity, current biomechanical models in the computation
community often assume only hyperelasticity. Fractional viscoelastic
constitutive models have been successfully used in literature to capture the
material response. However, the translation of these models into computational
platforms remains limited. Many experimentally derived viscoelastic
constitutive models are not suitable for three-dimensional simulations.
Furthermore, the use of fractional derivatives can be computationally
prohibitive, with a number of current numerical approximations having a
computational cost that is $ \mathcal{O} ( N_T^2) $ and a storage cost that is
$ \mathcal{O}(N_T) $ ($ N_T $ denotes the number of time steps). In this paper,
we present a novel numerical approximation to the Caputo derivative which
exploits a recurrence relation similar to those used to discretize classic
temporal derivatives, giving a computational cost that is $ \mathcal{O} (N) $
and a storage cost that is fixed over time. The approximation is optimized for
numerical applications, and the error estimate is presented to demonstrate
efficacy of the method. The method is shown to be unconditionally stable in the
linear viscoelastic case. It was then integrated into a computational
biomechanical framework, with several numerical examples verifying accuracy and
computational efficiency of the method, including in an analytic test, in an
analytic fractional differential equation, as well as in a computational
biomechanical model problem.
