Acta Polytechnica CTU Proceedings https://doi.org/10.14311/APP.2024.49.0013 Acta Polytechnica CTU Proceedings 49:13–19, 2024 © 2024 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague FRACTIONAL ORDER MODELS OF VISCOELASTIC POLYMERIC SOLIDS UNDERGOING LARGE DEFORMATIONS Barbora Hálkováa,∗, Michal Benešb a Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague, Czech Republic b Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mathematics, Thákurova 7, 166 29 Prague, Czech Republic ∗ corresponding author: barbora.halkova@fsv.cvut.cz Abstract. We present a fractional order model for nonlinear visco-hyperelastic solids taking into account large deformations. A three-field form of the Hu-Washizu principle is introduced to create a stable finite element method in the context of nearly incompressible dynamics. The β-method (a generalized midpoint rule) for time discretization is implemented into a variational finite element framework for efficient computing of numerical approximations to the initial boundary-value problem for hyperbolic equation of motion. Finally, a 2-D cantilever beam problem with a step end load is considered in order to demonstrate the algorithm. Keywords: Fractional viscoelasticity, large deformations, springpot, fractional calculus, finite element, integration algorithm. 1. Introduction Many polymers exhibit time dependent behavior some- where between purely elastic and purely viscous ma- terials. As a result, a large number of Kelvin-Voigt or Maxwell elements (and thus a large number of ma- terial parameters) are needed to be identified from experimental data to obtain a reasonably accurate description of mechanical response. On the other hand, a fractional calculus, i.e. the theory of deriva- tives and integrals of non-integer order, seems to be an efficient tool for the theoretical modelling of vis- coelastic materials [1]. Theoretical models based on the fractional calculus allows us to describe viscoelas- tic materials with significantly less parameters than the standard approach. For a deeper discussion on this issue we refer the reader to [2]. The fractional viscoelastic model at small strains was introduced e.g. in [3–7]. On the other hand, although rubbery polymers typically exhibit large deformations in en- gineering applications, much less attention has been given to fractional viscoelasticity in combination with the finite strain theory [8]. The present work pro- vides a computational framework for modelling the fractional viscoelastic behaviour of polymeric solids at finite strains in the context of nearly incompressible dynamics. Let the open set Ω0 be the reference configura- tion of a given (compressible or nearly incompress- ible) body at time t0. Here, Ω0 is described by a set of continuously distributed points X (parti- cles or material points) which occupy a region within the Euclidean space E3. In the absence of displace- ment discontinuities, a one-to-one deformation map ϕ : Ω0 × [0, T ] → E3 describing a motion exists, Figure 1. Motion of body Ω0. such that any displaced position at a current time t ∈ [0, T ], where [0, T ] ⊂ R+ denotes the time of interest, is determined as x = ϕ(X, t) with the differ- ence being the displacement field u(X, t) = x − X. Here, the position vector X = (X1, X2, X3) repre- sents the particle X in the reference configuration Ω0, X = X1E1 +X2E2 +X3E3, where (E1,E2,E3) de- fines an orthogonal base system with origin 0. Hence, we have x = ϕ(X, t) = X + u(X, t), see Figure 1. The deformation gradient F (X, t) is obtained as the gradient of this deformation map F = ∇0ϕ, where ∇0 is the gradient with respect to X, and the Jacobian of the deformation map is given by J(X, t) = det F . We also define the right Cauchy-Green deformation tensors C(X, t) = F T F . Let Ω0 be the reference configuration of the body of interest with the boundary ∂Ω0 and T > 0 be the time horizon. Let ∂Ωu, ∂ΩP be smooth open disjoint subsets of ∂Ω0 such that ∂Ω0 = ∂Ωu ∪ ∂ΩP , ∂Ωu ̸= ∅ 13 https://doi.org/10.14311/APP.2024.49.0013 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Barbora Hálková, Michal Beneš Acta Polytechnica CTU Proceedings and ∂Ωu ∩∂ΩP = ∅. Let f0 be the body force per unit of mass, a given vector field defined on Ω0 × (0, T ), u0 and v0 : Ω0 → R3 be the prescribed initial position and velocity. Further, let t̆0 be the given traction prescribed on ∂ΩP × (0, T ) and let the displacement ŭ0 be prescribed on ∂Ωu × (0, T ). The local form of the initial boundary value problem for the momentum equation is given by the following system [9]: ρ0 ∂2u ∂t2 = ∇0 · P + f0 in Ω0 × (0, T ), (1) u = ŭ0 in ∂Ωu × (0, T ), (2) P · N0 = t̆0 in ∂ΩP × (0, T ), (3) u(·, 0) = u0 in Ω0, (4) ∂u ∂t (·, 0) = v0 in Ω0, (5) where N0 is the field normal to ∂ΩP , P = P (X,F (X, t)) is the first Piola-Kirchhoff stress. 2. Constitutive relationships 2.1. Hyperelastic material In this contribution, we are interested to solve the sys- tem (1)–(5), in the context of nearly-incompressible material behavior (rubber like materials), which re- quires a special numerical treatment – like mixed methods. In particular, the deformation is split into a volumetric part represented by J and an isochoric part of the right Cauchy-Green deformation tensor Ĉ, Ĉ = J− 2 3 C. As a consequence, the split permits a different treatment of the incompressible part. The deformation gradient F together with its con- jugate first Piola-Kirchhoff stress measure P , will be retained in order to defined the basic material rela- tionship. The hyperelastic constitutive equation can be generally expressed as: P (F ) = ∂W (F ) ∂F . (6) Typically for mixed methods [10, 11], the stored energy function W is additively decomposed into distortional part Ŵ and dilatational part U , namely: W (F ) = Ŵ (C) + U(J). (7) Recall that C = F T F and J = det F . A tradi- tional nearly incompressible neo-Hookean potential Ŵ (C) = 1 2µ(trĈ − 3) and a simple volumetric func- tion U(J) = 1 2κ(J − 1)2 are used in this work, where µ is the shear modulus and κ denotes the bulk mod- ulus. These assumptions give us the possibility to decompose the stress tensor into pure shear and bulk responses. The calculation is straightforward, we get: P (F ) = 2F ∂Ŵ ∂C + 2F ∂U ∂J ∂J ∂C (8) Figure 2. The standard viscoelastic model. and ∂J ∂C = √ det C ∂C = 1 2JC−1. (9) To account properly for nearly incompressible ma- terial response and to avoid difficulties concerning “locking” of the finite element procedure, we employ a three-field de Veubeke-Hu-Washizu principle [10]. Additional variables entering the mixed three-field formulation represent a strain variable θ which is equivalent to J : θ = J, (10) and the hydrostatic pressure: p = ∂U ∂J ∣∣∣∣ J=θ . (11) In the weak formulation of our problem and the subse- quent finite element procedure, introduced hereafter, Equations (10) and (11) are satisfied in a weak sense. 2.2. Fractional viscoelasticity In a viscoelastic model the stress depends not solely on the current strain (elastic model), but it also depends on the entire strain history. The standard viscoelastic model consists ofN Maxwell chains coupled in parallel, see Figure 2. The present approach is based on the assumption that a viscous response is characterized by a set of rate constitutive equations, namely for a nonequilibrium stress Qk (in chain k) as an internal variable, k = 1, 2, . . . , N . The constitutive model can be written as a set of coupled equations [8, 9, 12]: S = 2∂W ∂C − N∑ k=1 Qk, P = F S, (12) ∂Qk ∂t + 1 τk Qk = 1 τk ( 2∂Ŵk ∂C ) , Qk(0) = 0, (13) where k = 1, 2, . . . , N, W = ∑N k=1 Wk, 14 vol. 49/2024 Fractional viscoelastic models at large deformations S represents the second (or symmetric) Piola- Kirchhoff stress tensor, Wk is the strain energy in chain k, τk is the relaxation time associated with each Maxwell chain. Classical theory of viscoelasticity employs the mod- els composed of rheological elements such as elastic springs and viscous dampers. Meanwhile, the frac- tional viscoelasticity introduces the springpot element together with the principles of fractional calculus. The fractional derivative of order α of a function u is de- fined as: Dαu(t) = 1 Γ(1 − α) d dt ∫ t 0 u(s) (t− s)α ds, (14) where 0 < α < 1, Γ is the gamma function. Replacing now the integer order derivative in (13) with a fractional order derivative we get: Dαk Qk + 1 ταk k Qk = 1 ταk k ( 2∂Ŵk ∂C ) , (15) Qk(0) = 0, (16) where k = 1, 2, . . . , N, ταk k can now be interpreted as the most probable relaxation time out of a continuous distribution of relaxation times. The fractional order of differentiation αk then plays the role of a distribution parameter for the corre- sponding distribution of relaxation times [13]. The fractional Maxwell chain is obtained as a parallel con- nection of N fractional Maxwell cells, see the scheme in Figure 3. 3. The weak formulation In the rest of the paper, just to simplify and shorten the presentation and avoid unnecessary technicali- ties, we will consider the case N = 1 (the fractional Maxwell cell). The simple model presented here can be straightforwardly extended to a setting with the par- allel connection of N cells (in the fractional Maxwell chain). Let W 1,2(Ω0) denote the Sobolev space of func- tions possessing square integrable derivatives and H(Ω0) := W 1,2(Ω0)3. Let us denote by St the dis- placement solution space at time t ∈ [0, T ] defined as: St = { u(·, t) ∈ H(Ω0) ∣∣∣∣ u(·, t) = ŭ0(·, t) on ∂Ωu } . Finally, we denote by Hu(Ω0) the linear space of ad- missible test functions or kinematically admissible Figure 3. The fractional Maxwell chain. variations, i.e., (virtual) displacements satisfying the homogeneous form of the essential boundary condi- tion (2) as: Hu(Ω0) = { v ∈ H(Ω0) ∣∣∣∣ v = 0 on ∂Ωu } . With these notations in hand, the weak form of our problem reads as follows: For sufficiently smooth data f0, u0, v0, t̆0 and ŭ0 find the displacement field u(·, t) ∈ St, pressure p(·, t) ∈ L2(Ω0), the volume filed θ(·, t) ∈ L2(Ω0) and the internal variable Q(·, t) ∈ L2(Ω0)3×3, such that u(·, 0) = u0, ∂u ∂t (·, 0) = v0 and Q(0) = 0 in Ω0 and: ∫ Ω0 ρ0 ∂2u ∂t2 · v dΩ0 + ∫ Ω0 Ŝ : [F T ∇0v] + pJC−1 : [F T ∇0v] dΩ0 − ∫ Ω0 Q : [F T ∇0v] dΩ0 = ∫ Ω0 f0 · v dΩ0 + ∫ ∂ΩP t̆0 · v dS0 (17) for all v ∈ Hu(Ω0) and all t ∈ [0, T ],∫ Ω0 [p− κ(θ − 1)]ψ dΩ0 = 0 (18) for all ψ ∈ L2(Ω0) and all t ∈ [0, T ] and:∫ Ω0 (J − θ)q dΩ0 = 0 (19) for all q ∈ L2(Ω0) and all t ∈ [0, T ]. Finally: DαQ + 1 τα Q = 1 τα Ŝ, τ > 0, α ∈ (0, 1], (20) holds almost everywhere in Ω0 and for all t ∈ [0, T ]. Here: Ŝ = 2∂Ŵ ∂C 15 Barbora Hálková, Michal Beneš Acta Polytechnica CTU Proceedings denotes the computational deviatoric second Piola- Kirchhoff stress. The Equation (17) denotes the weak form of the Equation of motion (1). The second Equation (18) yields the constitutive equation for the pressure p, see also Equation (11), and the third Equa- tion (19) reproduces the constraint condition (10). 4. Numerical algorithm Here we outline a general numerical solution scheme for the viscoelastic problem (extended to fractional viscoelasticity) within the context of the finite-element method. The point of departure in our developments is the weak formulation introduced in the preceding section. For simplicity, we will assume that the dis- placement boundary conditions are homogeneous, i.e.: u(·, t) = ŭ0(·, t) ≡ 0 on ∂Ωu. We begin by traditional discretization in space by the finite element method. Let us define: uh = ΨT u ũ, uh ∣∣∣∣ Ωe = NNu∑ i=1 Ψi u(x)ũi(t) (21) for all uh ∈ Hh u ⊂ Hu(Ω0), θh = ΨT θ θ̃, θh ∣∣∣∣ Ωe = NNθ∑ i=1 Ψi θ(x)θ̃i(t) (22) for all θh ∈ Lh ⊂ L2(Ω0), ph = ΨT p p̃, ph ∣∣∣∣ Ωe = NNp∑ i=1 Ψi p(x)p̃i(t), (23) for all ph ∈ Lh ⊂ L2(Ω0) and: Qh = ΨT QQ̃, Qh ∣∣∣∣ Ωe = NNQ∑ i=1 Ψi Q(x)Q̃i(t) (24) for all Qh ∈ [Lh]3×3 ⊂ L2(Ω0)3×3. Here we denoted by Hh u and Lh the finite element subspace of the space Hu(Ω0) and L2(Ω0), respectively. Let now 0 = t0 < t1 < · · · < tR = T be an equidis- tant partitioning of the time interval [0, T ] with the discrete time step ∆t, {tn}R n=0, ∆t = T R . For any function, vector-valued function or tensor ζ, we will use the approximation ζn ≈ ζ(tn) and introduce the notation: ζn+β = (1 − β)ζn + βζn+1, β ∈ [0, 1]. (25) Our goal is to develop time discretization schemes for which a discrete form of the problem (17)– (20) can be established. Replacing the second or- der derivative in (17) by the discrete derivative, ∂2u ∂t2 ≈ uh n+1−2uh n+uh n−1 ∆t , and incorporating the gen- eralized midpoint rule (25) into the weak formu- lation (17)–(19), the final result is the generalized β-scheme of the form: Ru = − ∫ Ω0 ρ0 ( uh n+1 − 2uh n + uh n−1 ∆t ) · δuh dΩ0 − β∆t ∫ Ω0 Ŝn−1+β : [ F T n−1+β(∇0 δu h) ] dΩ0 − β∆t ∫ Ω0 ph n−1+β Jn−1+βC−1 n−1+β : [ F T n−1+β(∇0 δu h) ] dΩ0 − β∆t ∫ Ω0 Qn−1+β : [ F T n−1+β(∇0 δu h) ] dΩ0 − (1 − β)∆t ∫ Ω0 Ŝn+β : [ F T n+β(∇0 δu h) ] dΩ0 − (1 − β)∆t ∫ Ω0 ph n+βJn+βC−1 n+β : [ F T n+β(∇0 δu h) ] dΩ0 − (1 − β)∆t ∫ Ω0 Qn+β : [ F T n+β(∇0 δu h) ] dΩ0 + β∆t ∫ Ω0 (f0)n−1+β · δuh dΩ0 + (1 − β)∆t ∫ Ω0 (f0)n+β · δuh dΩ0 + β∆t ∫ ∂ΩP (t̆0)n−1+β · δuh dS0 + (1 − β)∆t ∫ ∂ΩP (t̆0)n+β · δuh dS0 = 0 (26) for all δuh ∈ Hh u, Rθ = − ∆t ∫ Ω0 [ κβ ( θh n−1+β − 1 ) + κ(1 − β) ( θh n+β − 1 )] δθh dΩ0 √ + ∆t ∫ Ω0 [ βph n−1+β + (1 − β)ph n+β ] δθh dΩ0 = 0 (27) for all δθh ∈ Lh and: Rp = − ∆t ∫ Ω0 [βJn−1+β + (1 − β)Jn+β ] δph dΩ0 + ∆t ∫ Ω0 [ βθh n−1+β + (1 − β)θh n+β ] δph dΩ0 = 0 (28) for all δph ∈ Lh. Note that β serves as a parameter to control the implicitness of the algorithm. For β = 0 or β = 1 we get the explicit method. Setting β = 1 2 we obtain the implicit method. Next we present the algorithm for the integration of the rate Equation (20). First, the discrete frac- tional derivative of order α at time (n+ 1)∆t can be approximated by [14]: [DαQ]n+1 = 1 (∆t)α n∑ j=0 wj(α)Qn+1−j , (29) where the weights can be identified and calculated by the recursion relation below: w0(α) = 1, w1(α) = −α, · · · wj(α) = j − 1 − α j wj−1(α), · · · . 16 vol. 49/2024 Fractional viscoelastic models at large deformations Using (29) in the rate equation (20) we get: [DαQ]n+1 + 1 τα Qn+1 = 1 τα Ŝn. (30) Moreover, using (29) we can rewrite the approxima- tion of the fractional derivative as: [DαQ]n+1 = 1 (∆t)α Qn+1 + 1 (∆t)α n∑ j=1 wj(α)Qn+1−j . (31) Substituting (31) into (30) it is easy to see that Qn+1 can be computed as: Qn+1 =(∆t)α[w0(α)τα + (∆t)α]−1Ŝn − τα[w0(α)τα + (∆t)α]−1 n∑ j=1 wj(α)[Q]n+1−j . (32) The disadvantage of the fractional viscoelasticity is the nonlocal character of fractional derivatives. Nu- merical approximation requires the whole history of the internal variables in the preceding time steps to be saved and included in the calculation of the new time step. Letting α = 1 and the sum in Equation (32) simply becomes −Qn and the classical model follows. It is worth pointing out that the right hand side in (30) is taken from the preceding time step (a semi- implicit integrator). As a result, the discrete rate Equation (30) is decoupled from the system at the actual time step n+ 1 and the approximation of the internal variable Qn+1 can be computed first (at each discrete time step). Then, with Qn+1 in hand, the Newton-Raphson method is applied to solve the non- linear system (26)–(28). The derivation of a tangent stiffness tensor that is consistent with the integration procedure is essential to ensure a quadratic rate of convergence [15]. A consistent linearization for the set of non-linear equations given in Equations (26)–(28), about a configuration u, θ and p, is given by: Jg(y(i))∆y(i) = −g(y(i)), (33) y(i+1) = y(i) + ∆y(i), (34) where yT = (ũn+1, θ̃n+1, p̃n+1), gT = (Ru,Rθ,Rp), Jg =  Cuu 0 CT pu 0 Cθθ CT pθ Cpu Cpθ 0  , C□△ = ∂R□ ∂△i+1 . 5. Numerical example To illustrate the performance of the model, we con- sider a 2-D cantilever beam 1.0×0.4 m with a step end load, see Figure 4 and Figure 5. The beam is meshed uniformly with 4-node quadrilateral elements and par- titioned as shown. The parameters used for this prob- lem are: Young’s modulus E = 1.0 × 104 Pa, density Figure 4. 2-D cantilever beam. Figure 5. A step end load. Figure 6. Vertical displacement of the mid-point A on the free edge. Hyperelastic material was considered in this case (neglecting viscous effects, i.e. τ → +∞). ρ0 = 1.0 × 10−4 kg m−3, end load P0 = 20.0 N m−1 distributed over the free edge. The problem was integrated with a time step ∆t = 2.0 × 10−6 s. First, we considered the purely hyperelastic material ne- glecting viscous effects (τ → +∞). Figure 6 shows the motion of point A. Note that our results are in very close agreement with the results presented in [16]. In Figure 7 the vertical displacement of the mid- point A versus time is displayed for different values of α. As can be observed, the value of α clearly affects the results. Vertical displacements seem to decay faster with higher values of α. Figures 8 and 9 show the deformed cantilever for parameter α = 0.5. Due to the fact that the time step is limited by 17 Barbora Hálková, Michal Beneš Acta Polytechnica CTU Proceedings Figure 7. Vertical displacement of the point A. The influence of different values of α for fixed τ = 10−3 s. Figure 8. Motion of the beam at time t = 3.0 × 10−4 s. Figure 9. Motion of the beam at time t = 8.9 × 10−4 s. accuracy requirements rather than stability conditions, the explicit algorithm was used in our simulations (which is significantly faster as no iterations in each time step are needed). 6. Conclusion A fractional derivative visco-hyperelastic model for large and nearly incompressible deformations has been formulated based on irreversible thermodynamics with internal variables. The finite element framework is based on a three-field form of the Hu-Washizu prin- ciple to create a stable finite element method. The β-method (the generalized midpoint rule) for time discretization of the equation of motion and a specific semi-implicit approximation of fractional ODEs gov- erning the evolution of the internal variables enable us to partially decouple the elastic and viscous response to simplify and speed up the numerical algorithm. The consistent linearization of the resulting system of nonlinear equations is also briefly presented. The present work is our first step toward linking the fractional calculus and hyperelasticity. The dy- namic response of a 2-D cantilever beam is computed, including both geometrically and materially nonlinear effects. Acknowledgements This research was supported by SGS, project number SGS23/001/OHK1/1T/11, and the Czech Science Foun- dation, the grant No. 22-15553S. References [1] R. C. Koeller. Applications of fractional calculus to the theory of viscoelasticity. Journal of Applied Mechanics 51(2):299–307, 1984. https://doi.org/10.1115/1.3167616 [2] B. Hálková. Experimentální a numerické modelování PVB folie [In Czech; Experimental and numerical modelling of PVB foil]. Master’s thesis, Czech Technical University in Prague, 2024. [3] S. W. J. Welch, R. A. L. Rorrer, J. R. G. Duren. Application of time-based fractional calculus methods to viscoelastic creep and stress relaxation of materials. Mechanics of Time-Dependent Materials 3(3):279–303, 1999. https://doi.org/10.1023/A:1009834317545 [4] J. Padovan. Computational algorithms for FE formulations involving fractional operators. Computational Mechanics 2(4):271–287, 1987. https://doi.org/10.1007/BF00296422 [5] M. Enelund, L. Mähler, K. Runesson, B. L. Josefson. Formulation and integration of the standard linear viscoelastic solid with fractional order rate laws. International Journal of Solids and Structures 36(16):2417–2442, 1999. https://doi.org/10.1016/S0020-7683(98)00111-5 [6] S. Müller, M. Kästner, J. Brummund, V. Ulbricht. A nonlinear fractional viscoelastic material model for polymers. Computational Materials Science 50(10):2938–2949, 2011. https://doi.org/10.1016/j.commatsci.2011.05.011 [7] A. Schmidt, L. Gaul. Finite element formulation of viscoelastic constitutive equations using fractional time derivatives. Nonlinear Dynamics 29(1):37–55, 2002. https://doi.org/10.1023/A:1016552503411 [8] K. Adolfsson, M. Enelund. Fractional derivative viscoelasticity at large deformations. Nonlinear Dynamics 33(3):301–321, 2003. https://doi.org/10.1023/A:1026003130033 18 https://doi.org/10.1115/1.3167616 https://doi.org/10.1023/A:1009834317545 https://doi.org/10.1007/BF00296422 https://doi.org/10.1016/S0020-7683(98)00111-5 https://doi.org/10.1016/j.commatsci.2011.05.011 https://doi.org/10.1023/A:1016552503411 https://doi.org/10.1023/A:1026003130033 vol. 49/2024 Fractional viscoelastic models at large deformations [9] J. C. Simo, T. J. R. Hughes. Computational Inelasticity. Interdisciplinary Applied Mathematics Volume 7. Springer New York, USA, 1998. https://doi.org/10.1007/b98904 [10] J. C. Simo, R. L. Taylor, K. S. Pister. Variational and projection methods for the volume constraint in finite deformation elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 51(1–3):177–208, 1985. https://doi.org/10.1016/0045-7825(85)90033-7 [11] S. N. Atluri, E. Reissner. On the formulation of variational theorems involving volume constraints. Computational Mechanics 5(5):337–344, 1989. https://doi.org/10.1007/BF01047050 [12] J. C. Simo. On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects. Computer Methods in Applied Mechanics and Engineering 60(2):153–173, 1987. https://doi.org/10.1016/0045-7825(87)90107-1 [13] M. Enelund, G. A. Lesieutre. Time domain modeling of damping using anelastic displacement fields and fractional calculus. International Journal of Solids and Structures 36(29):4447–4472, 1999. https://doi.org/10.1016/S0020-7683(98)00194-2 [14] C. Lubich. Discretized fractional calculus. SIAM Journal on Mathematical Analysis 17(3):704–719, 1986. https://doi.org/10.1137/0517050 [15] J. C. Simo, R. L. Taylor. Consistent tangent operators for rate-independent elastoplasticity. Computer Methods in Applied Mechanics and Engineering 48(1):101–118, 1985. https://doi.org/10.1016/0045-7825(85)90070-2 [16] A. Prakash, K. D. Hjelmstad. A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics. International Journal for Numerical Methods in Engineering 61(13):2183–2204, 2004. https://doi.org/10.1002/nme.1136 19 https://doi.org/10.1007/b98904 https://doi.org/10.1016/0045-7825(85)90033-7 https://doi.org/10.1007/BF01047050 https://doi.org/10.1016/0045-7825(87)90107-1 https://doi.org/10.1016/S0020-7683(98)00194-2 https://doi.org/10.1137/0517050 https://doi.org/10.1016/0045-7825(85)90070-2 https://doi.org/10.1002/nme.1136 Acta Polytechnica CTU Proceedings 49:13–19, 2024 1 Introduction 2 Constitutive relationships 2.1 Hyperelastic material 2.2 Fractional viscoelasticity 3 The weak formulation 4 Numerical algorithm 5 Numerical example 6 Conclusion Acknowledgements References