% % internal variable discussion document % \documentclass[10pt, a4paper, fleqn]{article} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% \usepackage{/vol/PACK/tex_v3.14/lib/texmf/tex/tools/showkeys} %%% TeXdraw macros for diagrams... \usepackage{epsf} \usepackage{texdraw} \everytexdraw{\drawdim mm \linewd 0.2 \arrowheadtype t:F \arrowheadsize l:2 w:1} \def\thinlines{\linewd 0.008 } \def\normallines{\linewd 0.02 } \def\thicklines{\linewd 0.07 } % \def\VWall{ % vertical right facing wall: length = .3 \rmove(.3 -.3) \thicklines \rlvec(0 .6) \thinlines \rmove(-.1 -.6) \rlvec(.1 .1) \rmove(-.2 -.1) \rlvec(.2 .2) \rmove(-.2 -.1) \rlvec(.2 .2) \rmove(-.2 -.1) \rlvec(.2 .2) \rmove(-.2 -.1) \rlvec(.2 .2) \rmove(-.2 -.1) \rlvec(.2 .2) \rmove(-.2 -.1) \rlvec(.1 .1) \rmove(.1 -.3)\normallines } % \def\Spring#1{ % horizontal spring: length = 2.6 \rlvec(.3 0) \rmove(0 .4)\htext{\hbox to 2cm{\hfil#1\hfil}} \rmove(0 -.4) \rlvec(.1 .2) \rlvec(.2 -.4) \rlvec(.2 .4) \rlvec(.2 -.4) \rlvec(.2 .4) \rlvec(.2 -.4) \rlvec(.2 .4) \rlvec(.2 -.4) \rlvec(.2 .4) \rlvec(.2 -.4) \rlvec(.1 .2) \rlvec(.3 0) } % \def\DashPot#1{ % horizontal dashpot: length = 2.6 \rlvec(.81 0) \rmove(0 .7)\htext{\hbox to 1.3cm{\hfil#1\hfil}} \rmove(0 -.7) \rlvec(.6 0)\rmove(0 -.35)\rlvec(0 .7)\rmove(0 -.35) \rmove(-.4 -.4)\rlvec(.8 0)\rlvec(0 .8)\rlvec(-.8 0)\rmove(.8 -.4) \thinlines \rmove(-.3 -.3)\rlvec(0 .5)\rmove(.1 .1)\rlvec(0 -.5)\rmove(.1 -.1)\rlvec(0 .5) \rmove(.1 -.2)\normallines \rlvec(.81 0) } % \def\LabelledNode#1{ % node labelled vertically above \rlvec(.3 0)\fcir f:0 r:.08 \thinlines \rmove(0 .2) \rlvec(0 .6)\rmove(0 -.2) \htext{\hbox{\ #1}} \rmove(0 -.1)\ravec(1.2 0)\rmove(-1.2 -.5) \normallines \rlvec(.3 0) } % \def\LabelledRArrowNode#1{ % right arrow node labelled at right \ravec(.6 0)\htext{\hbox{\ #1}} } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % format control: \pagestyle{empty} \overfullrule=20pt \setlength{\parskip}{1.5ex}\setlength{\parindent}{2em} \setlength{\mathindent}{\the\parindent} \renewcommand{\baselinestretch}{1.05} \setlength{\textwidth}{13cm} \setlength{\textheight}{19.05cm} \setlength{\topmargin}{-1.3cm} \setlength{\headsep}{1cm} \setlength{\oddsidemargin}{1.45cm} \setlength{\evensidemargin}{0.6cm} %...etc % \newcommand{\xref}[1]{(\ref{#1})} \renewcommand{\emph}[1]{\textbf{#1}} \def\Apos{{\em A posteriori}\/} \def\apos{{\em a posteriori}\/} \def\Apri{{\em A priori}\/} \def\apri{{\em a priori}\/} \def\etal{{\em et al.}\/} \def\etc{{\em etc.}\/} \newfont{\Eusm}{eusm10 at 12pt}\newfont{\sEusm}{eusm10 at 8pt} \newfont{\Eufm}{eufm10 at 12pt} \newfont{\Amsa}{msam10 at 12pt} \newfont{\Amsb}{msbm10 at 12pt}\newfont{\sAmsb}{msbm10 at 8pt} \newcommand{\Real} {\mbox{\Amsb R}} \newcommand{\bold}[1]{\ensuremath{\mbox{\boldmath$#1$\unboldmath}}} \newcommand{\PDV}{\textsc{pdv}} \newcommand{\Insert}{\ \\\par\noindent Picture here\\\par\noindent} \def\Implies{\quad\Longrightarrow\quad} % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %%%%%% \title{\ \\\ \\A comparison of hereditary integral and internal variable approaches to numerical linear solid viscoelasticity % \footnote{This document is TR 97/2 at % \texttt{www.brunel.ac.uk/~icsrbicm/tech_rep.html}} } \author{Simon Shaw, M.~K.~Warby and J.~R.~Whiteman} \date{\textit{BICOM, Brunel University, Uxbridge, UB8 3PH, England.}} \maketitle \begin{quote} \textsc{abstract}\ \small This article first compares the mathematical models obtained when a problem of linear quasistatic solid viscoelasticity is modelled by: (i) a hereditary integral with Prony series; and, (ii) an evolution equation for internal strain variables. We then confirm that while the models are the same, the numerical approximations arising from them will, in general, be different. This suggests the need for \apos\ error estimates for the internal variable formulation so that viscoelasticity problems may be efficiently and reliably simulated. \end{quote} \begin{quote} \textbf{Acknowledgement}\ \small\itshape The authors are pleased to acknowledge the stimulation of their interest in internal variable methods in the viscoelasticity context derived from conversations with Dr A R Johnson. \end{quote} \thispagestyle{empty} \begin{center}\section{Introduction}\end{center} The purpose of this article is to discuss the similarities and differences between quasistatic viscoelasticity models based on constitutive relationships given in the form of: \begin{itemize} \item a hereditary integral; \item evolution equations for internal variables. \end{itemize} The hereditary integral form appears to be that most commonly employed throughout the viscoelasticity literature, see for example Refs \cite{JDF, GGa}, and has been used by the authors in, for example, Refs \cite{SWd, SWWb} to derive numerical schemes for the quasistatic linear viscoelastic problem. Recently however, A.~Johnson has suggested that the statement of the problem in terms of internal variables offers greater flexibility and practicality in implementing software for viscoelastic stress analysis. Accepting this to be the case it is then important to provide \apos\ error estimates for the resulting schemes in order to facilitate efficient and reliable computer simulation. For an introduction to A.~Johnson's work we refer for example to the article by Johnson and Tessler in Ref \cite{AJohnsonATesslerMaf}. In \cite{SWa, SWd} Shaw and Whiteman have employed the adaptive finite element techniques of K.~Eriksson \etal\ (e.g. Refs \cite{EJc,ErikEstepHansJohn}) and derived such error estimates for the Partial Differential Volterra (\PDV) Equations that model viscoelastic deformation when hereditary integrals are used to provide the constitutive relationship. The constitutive models based on hereditary integrals and internal variables are completely equivalent when the relaxation functions (see below) in the integrals are taken in the form of a \emph{Prony series}: % \begin{equation}\label{eq:prony1} \varphi(t) = \varphi_0 + \sum_{i=1}^N \varphi_i e^{-t/\tau_i}, \qquad\varphi_i > 0,\quad\tau_i\ge 0. \end{equation} % Initially we had hoped that this equivalence would mean that the numerical schemes based on these two alternative constitutive models could also be shown to be equivalent (in some sense). If this were the case then we would have aimed to show also that the \apos\ error estimates derived for the Volterra formulations could be carried over in a straightforward manner to A.~Johnson's internal variable formulation. It turns out that the numerical schemes are not equivalent and that they (probably) never can be. This is what we demonstrate below. The first section outlines the hereditary constitutive law, and the second then outlines the ``spring and dashpot'' model of viscoelasticity and illustrates how a constitutive equation may be derived in terms of an internal strain variable. It is a simple matter to show that this model is equivalent to the hereditary integral if the viscoelastic material is described by the Prony series \xref{eq:prony1}. We then outline discretization schemes for both of these models and show that they are, somewhat disappointingly, not the same. In this part we confine our attention to a simple model problem, although greater generality can be achieved if desired, see for example \cite{SWWb, AJohnsonATesslerMaf}. \begin{center}\section{Hereditary integral formulation}\end{center} In a time interval $\mathcal{J} := [0,T]$ suppose that the interior of a viscoelastic solid body occupies the region $\Omega\subset\Real^n$, and that it has a piecewise smooth boundary coinciding with $\partial\Omega$. The body is subjected to a system of body forces $\bold{f} = \{f_i(\bold{x},t)\}_{i=1}^n$ for $\bold{x}\in\Omega$ and $t\in\mathcal{J}$, and surface tractions $\bold{g} = \{g_i(\bold{x},t)\}_{i=1}^n$ for $\bold{x}\in\Gamma_N:=\partial\Omega\setminus\Gamma_D$ and $t\in\mathcal{J}$. Here $\Gamma_D$ is a time independent subset of $\partial\Omega$ with positive measure. We denote the displacement of the body from its equilibrium configuration by $\bold{u} = \{u_i(\bold{x},t)\}_{i=1}^n$. Assuming negligible inertia Newton's second law of motion gives for $t\in\mathcal{J}$ the following equilibrium equations with boundary data. \begin{eqnarray}\label{eq:eqeq} -\sum_{j=1}^n {\partial \sigma_{ij}\over \partial x_j} (\bold{x},t) & = & f_i(\bold{x},t),\qquad\mbox{ for }\bold{x}\in\Omega,\\ \bold{u}(\bold{x},t) & = & \bold{0}, \qquad\mbox{ for }\bold{x} \in\Gamma_D,\\ \widehat{n}_j(\bold{x})\sigma_{ij}(\bold{x},t) & = & g_i(\bold{x},t), \qquad\mbox{ for } \bold{x}\in\Gamma_N. \end{eqnarray} % In the above $\underline{\bold{\sigma}}=\{\sigma_{ij}\}_{i,j=1}^n$ is the (symmetric) stress tensor and $\widehat{\bold{n}} = \{n_j\}_{j=1}^n$ is a.e.~ the unit outward directed normal to $\Gamma_N$. In the linear theory the infinitesimal strain tensor $\underline{\bold{\varepsilon}}=\{\varepsilon_{ij}\}_{i,j=1}^n$ is defined in terms of the displacements as \begin{equation}\label{eq:strain} \varepsilon_{ij}(\bold{u}) := {1\over 2}\left( {\partial u_i\over\partial x_j} + {\partial u_j\over\partial x_i} \right). \end{equation} To close this problem we require a constitutive relationship between stress and strain. In linear elasticity theory this constitutive relationship is provided by Hooke's law: \[ \sigma_{ij} = D_{ijkl}\varepsilon_{kl}, \] where the fourth order tensor $\{D_{ijkl}\}$ describes the stiffness of the material. At this point it is convenient to alter our notation a little and think of $\underline{\bold{\sigma}}$ and $\underline{\bold{\varepsilon}}$ as vectors $\bold{\sigma}$ and $\bold{\varepsilon}$. Hooke's law may then be represented conveniently in matrix-vector form as \begin{equation}\label{eq:HookesLaw} \bold{\sigma} = \bold{D}\bold{\varepsilon}. \end{equation} The basic assumption made of linear viscoelastic materials is that the stress at the current time is a linear functional of not only the current strain but also of the entire strain history. One then introduces a time delay into the Hooke's tensor so that an increment in stress, $\Delta\bold{\sigma}(t)$, at the current time $t$ due to an increment in strain, $\Delta\bold{\varepsilon}(s)$, at time $s \le t$ is given by \begin{equation}\label{eq:delay} \Delta\bold{\sigma}(t) = \bold{D}(t-s)\Delta\bold{\varepsilon}(s). \end{equation} The assumption of linearity implies that the stress at current time is the initial stress $\bold{\sigma}(0)$ plus the sum of all stress increments at this time. Taking an appropriate limit of this superposition yields two equivalent forms of the linear constitutive relationship: \begin{eqnarray} \bold{\sigma}(t) & = & \bold{D}(0)\bold{\varepsilon}(t) - \int_0^t \bold{D}_s(t-s)\bold{\varepsilon}(s)\, ds, \label{eq:dispform}\\ \bold{\sigma}(t) & = & \bold{D}(t)\bold{\varepsilon}(0) + \int_0^t \bold{D}(t-s)\bold{\varepsilon}_s(s)\, ds. \label{eq:rateform} \end{eqnarray} In these the subscript $s$ denotes differentiation with respect to the ``historical time'' variable $s$ and (also in much of what follows) $x$ dependence has been suppressed. For the quasistatic problem it is most convenient to work with \xref{eq:dispform} since it generates ellipticity in the resulting \PDV\ equation. One now uses \xref{eq:dispform} with \xref{eq:eqeq} and \xref{eq:strain} to form the mathematical model as an elliptic equation with memory, or \emph{elliptic Volterra} equation. We refer to \cite{SWWDW} or \cite{SWWb} for details. \begin{center}\section{Internal variable formulation}\end{center} A perhaps more intuitive way of obtaining a constitutive relationship is to start with the physical observation that viscoelastic materials display the characteristics of both elastic solids and viscous fluids. The kinetics of these type of substances are modelled respectively by the spring and the dashpot. \begin{figure}[h] \caption{\scshape A Hookean (linear) spring: $\sigma = E\varepsilon$; $E$ is the spring stiffness} \bigskip\begin{center}\begin{texdraw} \drawdim cm \normallines \move(0 0) \VWall\rlvec(2 0)\Spring{$E$}\rlvec(2 0) \LabelledRArrowNode{$\varepsilon,\ \sigma = E\varepsilon$} \end{texdraw}\end{center} \end{figure} \begin{figure}[h] \caption{\scshape A Newtonian (linear) dashpot: $\sigma = \eta\displaystyle{d\varepsilon\over dt}$; $\eta$ is the viscosity} \bigskip\begin{center}\begin{texdraw} \drawdim cm \normallines \move(0 0) \VWall\rlvec(2 0)\DashPot{$\eta$}\rlvec(2 0) \LabelledRArrowNode{$\varepsilon,\ \sigma=\eta\displaystyle{d\varepsilon\over dt}$} \end{texdraw}\end{center} \end{figure} \noindent In these models the stress carried by the spring is proportional to the strain in the spring and is given by Hooke's law: $\sigma = E\varepsilon$. The stress carried in the dashpot is proportional to the strain rate and is given by Newton's law of viscosity: $\sigma = \eta\dot{\varepsilon}$. One then models a viscoelastic material by considering a notional system of springs and dashpots with independent stiffness and viscosity parameters. There are essentially two ways to connect a spring to a dashpot: in \emph{series} and in \emph{parallel}. These are the building blocks and are named the ``Maxwell'' and ``Voigt'' models. \section*{The Maxwell model} The Maxwell model is a series connection of a spring and dashpot. \begin{figure}[h] \caption{\scshape The Maxwell model} \bigskip\begin{center}\begin{texdraw} \drawdim cm \normallines \move(0 0) \VWall\rlvec(1 0) \Spring{$E$} \LabelledNode{$\varepsilon_S,\sigma_S$} \rlvec(1 0) \DashPot{$\eta$}\LabelledNode{$\varepsilon_D,\sigma_D$} \rlvec(1 0) \LabelledRArrowNode{$\varepsilon,\sigma$} \end{texdraw}\end{center} \end{figure} \noindent In this model $\varepsilon_S$ and $\sigma_S$ denote the strain and stress in the spring alone, and $\varepsilon_D$, $\sigma_D$ denote those in the dashpot alone. The total stress is given by $\sigma = \sigma_S = \sigma_D$ and the total strain by $\varepsilon = \varepsilon_S + \varepsilon_D$. Differentiating and using Hooke's and Newton's laws yield \[ {d\varepsilon\over d t} = {1\over E}{d \sigma_S\over d t} + {\sigma_D\over\eta} \Implies {d\sigma\over dt} + {\sigma\over\tau} = E{d\varepsilon\over dt}, \] where $\tau := \eta/E$ is the so-called \emph{relaxation time}. Using $\sigma(0) = E\varepsilon(0)$ this ODE is easily solved to give \[ \sigma(t) = Ee^{-t/\tau}\varepsilon(0) + E\int_0^te^{-(t-s)/\tau} \dot{\varepsilon}(s)\, ds, \] and this is essentially \xref{eq:rateform} with the scalar analogue of $\bold{D}$ given by the Prony series \xref{eq:prony1}. \section*{The Voigt model} Connecting the spring and dashpot in parallel yields the Voigt model. This time $\varepsilon_S=\varepsilon_D=\varepsilon$ and equilibrium demands that $\sigma = \sigma_S + \sigma_D$, hence % \begin{figure}[h] \caption{\scshape The Voigt model} \bigskip\begin{center}\begin{texdraw} \drawdim cm \normallines \move(0 0) \VWall\rlvec(1 0) \rmove(0 -1.5)\rlvec(0 3.0) \rlvec(.5 0) \DashPot{$\eta$} \LabelledNode{$\sigma_D,\varepsilon_D$} \rlvec(1.3 0) \move(1.3 -1.5) \rlvec(.5 0) \Spring{$E$} \LabelledNode{$\sigma_S, \varepsilon_S$} \rlvec(1.3 0) \rlvec(0 3.01)\rmove(0 -1.5)\rlvec(1 0)\LabelledRArrowNode{$\varepsilon,\sigma$} \end{texdraw}\end{center} \end{figure} % \[ \eta{d\varepsilon\over dt} + E\varepsilon = \sigma \Implies {d\varepsilon\over d t} + {\varepsilon\over \tau} = {\sigma\over \eta}. \] This gives the constitutive law in hereditary form as \[ \varepsilon(t) = e^{-t/\tau}\varepsilon(0) + {1\over\eta} \int_0^t e^{-(t-s)/\tau}\sigma(s)\, ds. \] \section*{The Maxwell solid} In his internal variable formulation A.~Johnson uses these basic building blocks in the Maxwell solid. \begin{figure}[h] \caption{\scshape The Maxwell solid} \bigskip\begin{center}\begin{texdraw} \drawdim cm \normallines \move(0 0) \VWall\rlvec(1 0) \rmove(0 -1.5)\rlvec(0 3.0)\rlvec(.5 0) \Spring{$E_1$}\LabelledNode{$\sigma^*,\varepsilon^*$} \rlvec(1 0)\DashPot{$\eta$} \LabelledNode{$\sigma_D,\varepsilon_D$} \move(1.3 -1.5)\rlvec(2.35 0) \Spring{$E_0$} \LabelledNode{$\sigma_S, \varepsilon_S$} \rlvec(2.35 0) \rlvec(0 3.01)\rmove(0 -1.5)\rlvec(1 0)\LabelledRArrowNode{$\varepsilon,\sigma$} \end{texdraw}\end{center} \end{figure} \noindent Here $E_0$ and $E_1$ are spring stiffnesses and $\sigma^*$, $\varepsilon^*$ are internal stress and strain variables. This time $\sigma^* = E_1\varepsilon^*$, $\varepsilon_D = \varepsilon-\varepsilon^*$ and $\sigma_S = E_0\varepsilon_S$. Also $\sigma^* = \sigma_D$ and this gives \[ E_1\varepsilon^* = \eta{d\over dt}(\varepsilon- \varepsilon^*) \Implies {d\varepsilon^*\over dt} + {\varepsilon^*\over\tau} = {d\varepsilon\over dt}, \] where now, and in all that remains, $\tau := \eta/E_1$. Solving this we get \begin{equation}\label{eq:ivhered} \varepsilon^*(t) = e^{-t/\tau}\varepsilon(0) + \int_0^t e^{-(t-s)/\tau}\varepsilon_s(s)\, ds. \end{equation} % Now, defining the \emph{stress relaxation function} \[ D(t) := E_0 + E_1 e^{-t/\tau} \] as the scalar analogue to the matrix $\bold{D}(t)$ in \xref{eq:dispform}, and using this in \xref{eq:ivhered} along with the relation \[ \sigma = \sigma_S + \sigma^* = E_0\varepsilon + E_1\varepsilon^* \qquad\mbox{(since $\varepsilon_S=\varepsilon$)}, \] gives \begin{eqnarray*} \sigma(t) & = & E_0\varepsilon(t) + E_1 e^{-t/\tau}\varepsilon(0) +\int_0^t E_1 e^{-(t-s)/\tau}\varepsilon_s(s)\, ds,\\ & = & D(0)\varepsilon(t) - \int_0^t D_s(t-s)\varepsilon(s)\, ds. \end{eqnarray*} Hence the internal variable and hereditary model (with Prony series) are completely equivalent. An important consequence of this equivalence is that the optimal long time stability estimates derived in \cite{SWWb} apply also to this internal variable formulation. Numerical schemes however will in general \emph{not} be equivalent, we demonstrate this below with the aid of a model problem. Note that according to A.~Johnson, \cite{AJohnsonPersonal}, it is essential that the spring be placed before the dashpot in this model. Reversing the order results in the differential equation \[ \eta{d\varepsilon^*\over dt} + E_1\varepsilon^* = E_1\varepsilon, \] which (apparently) is not suited to an internal variable formulation. \begin{center}\section{Numerical schemes}\end{center} Consider the one space-dimensional problem: \begin{equation}\label{eq:1dmodel} -\sigma_x(x,t) = f(x,t)\qquad\mbox{ in }\Omega := (0,1), \end{equation} and for $t\in\mathcal{J}$, with the boundary data \[ u(0,t) = 0\qquad\mbox{ and }\qquad\sigma(1,t) = g(t). \] We use the constitutive model ($x$ dependence suppressed) \begin{eqnarray} \sigma(t) & = & E_0\varepsilon(t) + E_1\varepsilon^*(t),\label{eq:1dinternal}\\ & = & D(0)\varepsilon(t) - \int_0^t D_s(t-s)\varepsilon(s)\, ds, \label{eq:1dhistory} \end{eqnarray} where the internal strain variable satisfies \begin{equation}\label{eq:1dODE} {\partial\varepsilon^*\over\partial t} + {\varepsilon^*\over\tau} = {\partial\varepsilon\over\partial t}, \end{equation} and we take $u^*(x,0) = u(x,0)$ where the internal displacement variable $u^*$ is defined in a natural way as $u_x^* := \varepsilon^*$. We now use \xref{eq:1dinternal} in \xref{eq:1dmodel} to give this problem a weak formulation: \\ for each $t$, find $u\in V:=\{v\in H^1(0,1):v(0)=0\}$ such that \[ a(u(t),v) = F(t,v) - \Phi(u^*(t), v)\qquad\forall v\in V. \] Here: \begin{eqnarray*} a(w,v) & := & \int_0^1 E_0 w_x v_x\, ds,\\ F(t,v) & := & \int_0^1 f(t) v\, dx + g(t) v(1),\\ \Phi(w,v) & := & \int_0^1 E_1 w_x v_x \, dx. \end{eqnarray*} By replacing $V$ with a finite element subspace $V^h$ we now pose the semidiscrete problem: \\ find $U:\mathcal{J}\rightarrow V^h$ such that \[ a(U(t), v) = F(t,v) - \Phi(U^*(t), v), \] where $U^*$ satisfies $U^*(0) = U(0)$ and \[ {d U^*\over dt} + {U^*\over\tau} = {d U\over dt} \qquad\mbox{ at every } x\in\Omega. \] This problem can of course be written in matrix form: \begin{equation}\label{eq:system} \bold{A}\bold{U}(t) = \bold{F}(t) - \bold{A}^*\bold{U}^*(t) \qquad\mbox{ with }\qquad {d \bold{U}^*\over dt} + {\bold{U}^*\over\tau} = {d \bold{U}\over dt}, \end{equation} and $\bold{U}^*(0) = \bold{U}(0)$. \begin{center}\section{The internal variable model}\end{center} A simple discretization of the system \xref{eq:system} consists of applying the backward Euler (discontinuous Galerkin) method to the ODEs. Here $\bold{U}$ and $\bold{U}^*$ are assumed to be piecewise constant on each time interval $\mathcal{J}_i := (t_{i-1},t_i)$, and given by $\bold{U}_i\approx \bold{U}\big\vert_{\mathcal{J}_i}$ and $\bold{U}^*_i\approx \bold{U}^*\big\vert_{\mathcal{J}_i}$. The ODE set then becomes \[ \bold{U}_i^* = \beta\bold{U}_{i-1}^* + \beta(\bold{U}_i - \bold{U}_{i-1}) \qquad\mbox{ where }\qquad \beta := \left(1+{k\over\tau}\right)^{-1}, \] and $k = t_i-t_{i-1}$ is the (constant) time step. Using $\bold{U}^*_0 = \bold{U}_0$ we can solve this and get \[ \bold{U}^*_i = \beta^i(\bold{U}_1-\bold{U}_0) + \beta^{i-1}(\bold{U}_2-\bold{U}_1) + \cdots + \beta(\bold{U}_i-\bold{U}_{i-1}) + \beta^i\bold{U}_0, \] for $i=1,2,3,\ldots$. Notice that for each $t_i$ \[ \lim_{k\rightarrow0} \beta^i = e^{-t_i/\tau}. \] Using this representation of $\bold{U}_i^*$ in \xref{eq:system} gives \[ \bold{A}\bold{U}_i = \bold{F}(t_i) - \bold{A}^*\left( \beta^i(\bold{U}_1-\bold{U}_0) + \beta^{i-1}(\bold{U}_2-\bold{U}_1) + \cdots + \beta(\bold{U}_i-\bold{U}_{i-1}) \right), \] which after some rearrangement yields, \begin{equation}\label{eq:compiv} (\bold{A}+\beta\bold{A}^*)\bold{U}_i = \bold{F}(t_i) - (\beta-1)\bold{A}^*\left( \beta\bold{U}_{i-1} + \beta^2\bold{U}_{i-2} + \cdots + \beta^{i-1}\bold{U}_{1} \right). \end{equation} Now let's turn our attention to the hereditary integral model. \begin{center}\section{The hereditary model}\end{center} Using the hereditary constitutive relationship \xref{eq:1dhistory} in \xref{eq:system} results in the system of second-kind Volterra equations. \[ (\bold{A}+\bold{A}^*)\bold{U}(t) = \bold{F}(t) + {1\over\tau}\int_0^t e^{-(t-s)/\tau}\bold{A}^*\bold{U}(s)\, ds. \] To effect a comparison between the numerical schemes we again approximate $\bold{U}$ by a piecewise constant function. This gives \begin{eqnarray*} (\bold{A}+\bold{A}^*)\bold{U}_i & = & \bold{F}(t_i) + \sum_{j=1}^i\bold{A}^*\bold{U}_j\ \int_{t_{j-1}}^{t_j} {e^{-(t_i-s)/\tau}\over\tau}\, ds,\\ & = & \bold{F}(t_i) + (1-\gamma)\bold{A}^*\bold{U}_i +\bold{A}^*\sum_{j=1}^{i-1}(1-\gamma)\gamma^{i-j}\bold{U}_j, \end{eqnarray*} where \[ \gamma := e^{-k/\tau}\Implies \gamma^i = e^{-t_i/\tau} \quad\mbox{ and }\quad\gamma^{i-j} = e^{-(t_i-t_j)/\tau}. \] Now solving for $\bold{U}_i$ and using \[ \bold{A} + \bold{A}^* - (1-\gamma)\bold{A}^* = \bold{A} + \gamma\bold{A}^*, \] gives \begin{equation}\label{eq:compint} (\bold{A}+\gamma\bold{A}^*)\bold{U}_i = \bold{F}(t_i) -(\gamma-1)\bold{A}^* \left( \gamma\bold{U}_{i-1}+\gamma^2\bold{U}_{i-2} + \cdots + \gamma^{i-1}\bold{U}_1 \right). \end{equation} This will be identical to \xref{eq:compiv} if and only if $\gamma^i = \beta^i$ for $i=1,2,3,\ldots$. This is only true when $k\rightarrow 0$. Since the ODE solver used in the internal variable formulation will never sample values of the exponential functions in the Prony series, while any quadrature rule used in the hereditary formulation in general must, we conclude that numerical approximations to these models will never be equivalent. Since A.~Johnson in \cite{AJohnsonPersonal} claims the internal variable formulation to be superior in terms of flexibility of implementation, our conclusions suggest that the development of adaptive solution schemes predicated on \apos\ error estimates should be a high priority to practitioners working in this field. \begin{center}\section{Conclusion}\end{center} In closing we wish to suggest an appropriate framework within which one could set up and analyze the internal variable model. If we identify stresses due only to elastic effects as $\underline{\bold{\sigma}}^E$, and stresses due to viscoelastic effects only as $\underline{\bold{\sigma}}^*$ then we can write \[ \sigma_{ij}^E = D_{ijkl}\varepsilon_{kl} \qquad\mbox{ and }\qquad \sigma_{ij}^* = D_{ijkl}^*\varepsilon_{kl}^*, \] where $\underline{\bold{\varepsilon}}$ is the tensor of actual strains present in the viscoelastic body, and $\underline{\bold{\varepsilon}}^*$ is the tensor of internal strains. Also these strains, \[ \varepsilon_{ij} := {1\over 2}\left( {\partial u_i\over\partial x_j} + {\partial u_j\over\partial x_i} \right) \qquad\mbox{ and }\qquad \varepsilon_{ij}^* := {1\over 2}\left( {\partial u_i^*\over\partial x_j} + {\partial u_j^*\over\partial x_i} \right), \] satisfy the differential constitutive equations: \[ {\partial\varepsilon_{ij}^*\over\partial t} +{\varepsilon_{ij}^*\over \tau} = {\partial \varepsilon_{ij}\over\partial t}. \] Newton's second law then gives the equations of motion as \[ \varrho {\partial^2 u_i\over\partial t^2}(t) - \sum_{j=1}^n {\partial\sigma_{ij}^E\over\partial x_j}(t) = \tilde{f}_i(t) \quad := \quad f(t) + \sum_{j=1}^n {\partial\sigma_{ij}^*\over\partial x_j}(t), \] where $\varrho$ is the mass density. As noted by A.~Johnson, the advantage here is that by replacing the actual load $\bold{f}(t)$ by $\tilde{\bold{f}}(t)$, which is modified by internal viscoelastic effects, standard \emph{elasticity software} can be applied to either the resulting dynamic or static problem. Indeed we might cautiously speculate that when large strains are present in the problem one might even be able to ``blend'' geometrically nonlinear elasticity models with linear (or even nonlinear) constitutive models. If so, then a similar blending might also be applied to error analysis and yield a very desirable framework in which to study these kind of problems. However, we should emphasize that these remarks are only \emph{very cautious speculation}. \bibliographystyle{plain} \bibliography{research} %%%%%% \end{document}