Calculating time derivatives of stress ($σ$) in the constitutive equation of a viscoelastic model without explicitly remembering the prior stress? I am trying to understand how to work with viscoelastic models. There are numerous simple viscoelastic models (described at that link). In each, there is a "constitutive relationship" between stress ($σ$) and strain ($Ɛ$).
One of the simpler examples to illustrate is the Kelvin-Voigt model. Here the constitutive relationship is easy to work with for my purposes because stress is only expressed in terms of strain which can be calculated geometrically at any point in time:
$$σ = EƐ + η\dot{Ɛ}$$
But when you get to higher complexity models like for example the "Standard Linear Model," it becomes less simple as current stress depends on the rate of change of stress:
$$σ = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ}$$
As you can see current stress depends on the prior stress as it requires the time derivative of stress to solve the current stress.
I understand there is a way to express the time derivatives like $\dot{σ}$ and $\ddot{σ}$ in terms of derivatives of strain instead. This would solve my problem with how I need to use these equations.
I don't understand how to do this however. I see someone did this in the first comment of the thread replying to this question, but I don't understand how it was done. They posted that for the three element diagram above, stress could be expressed both as:
$$σ = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ}$$
$$σ = E_1ϵ  + \frac{E_2ηϵ\dot{ϵ}}{E_2ϵ + η\dot{ϵ}}$$
I am told the second expression is actually wrong. However, I still believe there must be some way to remove the derivatives of stress in the solution for stress and express it all in terms of strain/geometry.
 A: Diagrams like these are used to build complex constitutive models from simple elements, very analogously to electric circuits built from resistors and capacitors. Parallel elemnts have the same deformation, and theri stresses add up; for series elements, it's the other way around. (Actually, there is an exact equivalence if you interchange stress/voltage, deformation rate/current, spring/capacitor, dashpot/resistor and parallel/series placement.)
So for example, the inline combination of spring and dashpot (the bottom leg of your second diagram) has a spring with relation ($\sigma_{E_2}=E\epsilon_{E_2}$) and a dashpot ($\sigma_\eta=\eta\dot\epsilon_\eta$). Since they're inline, the have the same stress ($\sigma_{E_2}=\sigma_\eta=\sigma_\text{bottom}$), and the deformations add up, $\epsilon_{E_2}+\epsilon_\eta=\epsilon_\text{bottom}$. Taking the time derivative of the last equation and inserting the previous two, you get $$\dot\epsilon_\text{bottom}=\dot\epsilon_{E_2}+\dot\epsilon_\eta=\frac{1}{E_2}\dot\sigma_\text{bottom} + \frac{1}{\eta}\sigma_\text{bottom}\,,\tag{*}$$
and voilá, you have one complete equation for the spring-dashpot inline combination. (Compare that to the parallel case (your first eqution), where you have $\epsilon$ and $\dot\epsilon$, but only $\sigma$.)
Now for the full three-element model: This is a parallel combination of the top and bottom parts, so $\epsilon_\text{top}=\epsilon_\text{bottom}=\epsilon$ and $\sigma_\text{top}+\sigma_\text{bottom}=\sigma$. For the bottom half, we have justr derived the relation $(*)$; for the top, we have $\sigma_\text{top}=E_1\epsilon$.
Now, we just need to eliminate $\sigma_\text{top/bottom}$ in favour of $\sigma$. But that's straightforward: In $(*)$, express $\sigma_\text{bottom}=\sigma-\sigma_\text{top}$, and use $\sigma_\text{top}=E_1\epsilon$. After some rearrangement, you get
$$\left(\frac{E_1+E_2}{E_2}\right)\dot\epsilon +\frac{E_1}{\eta}\epsilon= \frac{1}{E_2}\dot\sigma + \frac{1}{\eta}\sigma\,,\tag{**}$$
which is equivalent to your expression.
With the same general procedure, you can derive expressions for other combinations. Of course, so far, everything is linear. This implies that if you have two solutions for $\sigma$ and $\epsilon$, you can add them and get a new solution. Note that the expression by @Chemomechanics in the comments above is not linera, so it cannot be corretc for this type of circuit. (The expression gives the stress as a function of the strain, $\sigma=\sigma(\epsilon)$, but a sum of strains does not give the sum of stresses, $\sigma(\epsilon_1+\epsilon_2)\neq\sigma(\epsilon_1)+\sigma(\epsilon_2)$.)
In general, for circuits with dashpots, you cannot get rid of the time derivatives because the history is important, so you have to specify the history of stess or strain and then solve for the other.
There are more possible elements, such as

*

*a slider (St. Venant) element, which corresponds to ideal plastic deformation (no motion for stresses below a threshold, arbitrary movement at the threshold, or

*a mass, which wants to keep moving (dual to the inductor in the mechanic/electric duality), which is often neglected for continuum mechanics when we're dealing with slow processes.

On the other hand, you could add nonlinear generalisations, i.e. dashpots with a more complicated relation of strain rate and stress, to better describe real materials.
As a final note to answer the recenntly-added part of the question:
The constititutive relation $(**)$ is used together with the equation of motion (EOM), which in the "supporting article" is Eq. (2). The EOM contains $\sigma$ and a spatial derivative $\sigma{,x}$, together with terms involving the displacement $w$ and various derivatives, in this case up to second order in space and time,
$$ (\text{stuff})\cdot\sigma+ (\text{more stuff})\cdot\sigma_{,x}= (\text{stuff with $w$ and derivatives up to second order}) \tag{#}$$
To try and combine the EOM and Eq. $(**)$, you can rewrite Eq. $(**)$ as
$$\underbrace{\left[\left(\frac{E_1+E_2}{E_2}\right)\frac{\text d}{\text d t} +\frac{E_1}{\eta}\right]}_{\displaystyle \Xi}\epsilon= \underbrace{\left[\frac{1}{E_2}\frac{\text d}{\text d t} + \frac{1}{\eta}\right]}_{\displaystyle \Gamma}\sigma\,.\tag{##}$$
This is just rewriting the derivatives with operators and rearranging the terms into two larger operators. Now you can apply the operator $\gamma$ to Eq. $(\#)$, i.e. multiply both sides with $1/\eta$,  and add the time derivative times $(E_1+E_2)/E_2$. Then you will have a term $\Gamma \sigma$, i.e. the rhs of Eq. $(\#\#)$, which you can turn into an expression in terms of $\epsilon$. It's not clear to me right now how you would not end up up with terms where the derivatives act on the terms I've calles $(\text{stuff})$ in Eq. $(\#)$, but maybe you can rearrange and get rid of them -- the authors apparenty find some combinations involving the speed $c$. In any case, it makes sense that Eq. (14) in the paper now contains derivatives up to fourth order.
A: A linear general viscoelastic model has the form:
$$A^{(0)}\sigma+A^{(1)}\frac{\partial \sigma}{\partial t} + \dots +
A^{(n_1)}\frac{\partial^{n_1} \sigma}{\partial t^{n_1}} =
B^{(0)}\varepsilon + B^{(1)}\frac{\partial \varepsilon}{\partial t} + \dots +
B^{(n_2)}\frac{\partial^{n_2} \varepsilon}{\partial t^{n_2}}$$
Usign Laplace transform, you get:
$$\mathcal{L}[\sigma]\left(\sum_{k=0}^{n_1}A^{(k)}s^k\right)  - \left(\sum_{k=0}^{n_1}\sum_{j=1}^{k} A^{(k)}\frac{d^{(j-1)}\sigma(0)}{dt^{(j-1)}}s^{k-j}\right) =\\ \mathcal{L}[\varepsilon]\left(\sum_{k=0}^{n_2}B^{(k)}s^k\right)  - \left(\sum_{k=0}^{n_2}\sum_{j=1}^{k} B^{(k)}\frac{d^{(j-1)}\sigma(0)}{dt^{(j-1)}}s^{k-j}\right) $$
This allows us to write the relationship between the Laplace transform of stress and strain as:
$$A(s)\mathcal{L}[\sigma] = B(s)\mathcal{L}[\varepsilon] + C(s)$$
This allows you to write as you wish:
$$\sigma(t) = E\varepsilon(t) + \int_0^t F(t-\tau)\dot{\varepsilon}(\tau)d\tau, \\ \varepsilon(t) = \frac{\sigma(t)}{E} + \int_0^t K(t-\tau)\dot{\sigma}(\tau)d\tau$$
