Finite Element Formulation of Fractional Constitutive Laws Using the Reformulated Inﬁnite State Representation

: In this paper, we introduce a formulation of fractional constitutive equations for ﬁnite element analysis using the reformulated inﬁnite state representation of fractional derivatives. Thereby, the fractional constitutive law is approximated by a high-dimensional set of ordinary differential and algebraic equations describing the relation of internal and external system states. The method is deduced for a three-dimensional linear viscoelastic continuum, for which the hydrostatic and deviatoric stress-strain relations are represented by a fractional Zener model. One- and two-dimensional ﬁnite elements are considered as benchmark problems with known closed form solutions in order to evaluate the performance of the scheme.


Introduction
The mechanical behavior of viscoelastic materials is characterized by a combined viscous and elastic material response to external loads and displacements depending on time or frequency. It is often described by complex rheological models consisting of spring-dashpot networks, which usually lead to a large number of parameters and poor extrapolation properties over several decades in time or frequency. It has been shown, though, that additional fractional-order elements, so called springpots, may lead to improved approximations, even with few parameters [1,2]. These elements are related to integrals and derivatives of non-integer order, treated within the theory of fractional calculus [3][4][5][6].
Fractional derivatives and integrals arise in viscoelasticity, when their integral kernels are used to model long-term creep and relaxation processes in the time domain. Early contributions identify a power-law behavior of viscoelastic media, see [7][8][9][10][11]. Later, this empirical observation has been connected directly to the theory of fractional calculus, see [1,12,13] and the physical consistency of fractional viscoelastic models has been discussed in [14,15]. In most applications of fractional viscoelasticity, the mechanical behavior of rubber and polymers is investigated, see e.g., [2]. Recently, fractional viscoelastic models have been used to model asphalt and concrete creep, see [16][17][18], respectively.
A key problem of fractional calculus in applications is the efficient numerical solution of fractional-order differential equations. Unfortunately, the non-locality of fractional operators usually leads to computationally expensive implementations. One possible way to overcome these difficulties is the approximation of fractional-order differential equations by high-dimensional ordinary differential equations. This approximation is based on the diffusive or infinite state representation of fractional derivatives as introduced in [19,20] and elaborated in [21][22][23][24]. The associated numerical schemes (in the following called infinite state schemes) do not require to process the entire solution history in every time step, which leads to reduced computational costs and memory requirements, whereas the accuracy depends on the chosen quadrature and order of differentiation [25][26][27][28][29].
In order to simulate the mechanical behavior of arbitrary viscoelastic structures described by fractional constitutive laws, an associated finite element formulation is required, which has to be combined with a suitable solver for fractional-order differential equations. A first formulation for linear fractional viscoelastic systems is given in [30] by using the Laplace domain. This approach has been generalized to the time domain in [31,32]. Hereto, the equations of motion are transformed to order 2 + α and solved in terms of a fractionalorder state space vector. The authors in [33] derive the three-dimensional generalization of a fractional viscoelastic constitutive law for the isotropic case and introduce a finite element approach solving fractional-order equations of order 2 using the Grünwald-Letnikov scheme. The formulation is generalized for the anisotropic case in [34]. Both methods use an internal variable approach, given in general form in [35]. Another formulation, see [36], which also uses the Grünwald-Letnikov scheme, incorporates the constitutive law in a slightly different way. The method is adapted and improved in [2]. All authors mentioned so far consider the geometrically linear case. Formulations in the context of large deformations are given in [37] and, recently in [38,39]. Both of the latter contributions use an infinite state scheme similar as in [40] for the time integration. An additional accuracy optimization is included in [39].
A novel improved infinite state scheme was presented in [41], but is not yet accessible for a finite element formulation of fractional constitutive laws. It is the aim of this paper to fill that gap and combine the scheme in [41] with the finite element method (FEM). The method is based on the reformulated infinite state representation, which leads to a non-singular integration kernel that asymptotically decays to zero and hence results in small errors in an associated quadrature, called reformulated infinite state scheme (RISS). The finite element formulation is fully derived for an isotropic three-dimensional continuum described by a fractional Zener model (distinguishing between hydrostatic and deviatoric stress and strain components), being a simple but archetype fractional rheological model. The method is tested on a one-and a two-dimensional example with closed form exact solutions, allowing for a proper error analysis. Moreover, the characteristic structure of the reformulated infinite state representation is interpreted in terms of viscoelastic quantities of a springpot.
The necessary theoretical preliminaries of fractional calculus are given in Section 2, including the (reformulated) infinite state representation (Sections 2.2 and 2.3) and the scheme RISS (Section 2.4). The formulation of fractional viscoelasticity is provided in Section 3. Particularly, in Section 3.1, the terms and quantities of classical linear viscoelasticity are introduced and connected in Section 3.2 to the fractional-order case, especially the springpot and the fractional Zener model. The main results regarding FEM can be found in Section 4. It contains a three-dimensional formulation of the fractional Zener model (Section 4.1) and its implementation in the finite element method (Section 4.2).
The numerical examples are studied in Section 4.3. All results are discussed in Section 5.

Definitions and Properties
The following concepts and properties are well-known results from fractional calculus, see e.g., [3][4][5][6]. Let T > 0 and x ∈ L 1 (−∞, T] be an integrable function. The Liouville-Weyl fractional integral of x of order α ∈ (0, 1) is introduced as [42] where Γ represents the Euler Gamma function It can easily be checked that (1) simplifies to a plain primitive of x at the limit α → 1.
To examine the limit for α → 0, consider the following. If, additionally, the function x is absolutely continuous (in the sense of [6] (Chapter 2, Definition 6.1)) and x andẋ vanish in an appropriate manner in the negative time limit, we obtain by partial integration of (1) The representation (3) shows that (1) is a plain identity for α → 0. Both limiting cases together provide the intuition that the fractional integral for α ∈ (0, 1) is something in between the identity and the first-order primitive. Furthermore, the Liouville-Weyl fractional derivative of Caputo type of order α ∈ (0, 1) of an absolutely continuous function x is given as [42] which is (one possible definition for) a left inverse of the operator I α , i.e., see [6] (Chapter 2, § 6.2). Similar to the fractional integral, the fractional derivative interpolates between identity and first-order derivative. To see this, consider fractional derivatives of the function x(t) = t χ (0,∞) (t) in Figure 1, where χ (0,∞) is the unit function vanishing on (−∞, 0]. Further details regarding fractional integrals and derivatives are given, e.g., in [3][4][5] and especially for functions with unbounded domains in [6]. In order to use fractional calculus in viscoelasticity, certain classes of fractional-order differential equations have to be solved. For the applications in this article, we investigate the type of problem where Θ(t) is the Heaviside step function and γ, δ, x 0 ∈ R, γ = 0. It represents a linear functional equation with given initial function ϕ. By introducing the fractional Caputo derivative with zero initialization time and using the relationship the problem (6) can be rewritten as the fractional initial value problem The solution of (8) contains the so-called one-parameter Mittag-Leffler function see [43] and is through [3] (Theorem 7.2) given by
such that Accordingly, an infinite state representation of the fractional derivative C D α for α ∈ (0, 1) is given by with infinite states z(η, t). The two kinds of infinite states are related by the following equations, which may be obtained by using the solutions of the differential Equations (11b) and (17b) with initial values (11c), (17c) and partial integration as whenever x is bounded. The infinite state representation translates fractional integrals and derivatives to integer-order at the cost of a continuum of state variables. Correspondingly, the history of the function x is transferred to initial conditions of the infinite states (11c) and (17c). The approximation of the improper integrals of infinite states (11a) and (17a) by sums transforms fractional-order differential equations to ordinary differential equations or differential algebraic equations and indicates a possible starting point for the formulation of numerical schemes to solve fractional-order problems, see [41] (Section 3.2).

Reformulated Infinite State Representation
The reformulated infinite state representation is introduced as a rearrangement of (17a) aiming for an improved numerical integration of infinite states. To obtain the reformulation, the integrand in (17a) is expanded by the term λ 2 + ω 2 for a fixed ω > 0, which leads together with (11b), (17b) and (18) to The terms in (19) can be simplified by using the following properties proven in [41,45,46]. hold.
Consequently, (19) can be expressed as The advantage of the reformulated infinite state representation (22) is the new kernel with parameter ω > 0 which is integrable in (0, ∞) and fulfills Hence, the improper integrals in (22) tend to be more accurately computable by quadrature rules in comparison to the formulation (17a). Moreover, (22) contains only first-order derivatives of x and the infinite states Z(η, ·) and z(η, ·), being key to the subsequent numerical scheme, which is based on the solution of high-dimensional ODEs. For a fixed value of α, the function K ω (α, ·) has a maximum at Hence, the position of η max may be adjusted by the magnitude of ω. In the following, the abbreviation K(α, ·) := K ω (α, ·) is used, where ω fulfills

Reformulated Infinite State Scheme
The key idea of the reformulated infinite state scheme (RISS) introduced in [41] is to approximate the integrals of the infinite states Z andż in (22) by sums of a finite number of states performing a composite Gaussian quadrature We choose η 0 = 0 and η 1 , . . . , η K logarithmically equidistant in (0, ∞) and perform a Gaussian quadrature in each interval (η 0 , η 1 ), . . . , (η K−1 , η K ). Thereby, the shifted Gauss-Legendre nodes and weights are related to the standard Gauss-Legendre nodes and weights To abbreviate, denote the infinite states as Accordingly, in a differential equation (such as (6)), a fractional derivative C D α x(t), α ∈ (0, 1) can be approximated with (22) and (28) as together withŻ related to (11b) and (17b). Hence, the originally fractional-order problem may be approximated by an ordinary differential equation, which together with (34) can be solved by a standard ODE solver. Thereby, the initial function of the original problem has to be translated to initial values Z k,j (0), z k,j (0) for the approximating ODE through A detailed description and further evaluation of the scheme RISS is given in [41].

Introduction to Linear Viscoelasticity
The mechanical behavior of a deformable body is, aside from the general principles of mechanics (i.e., equilibrium and compatibility equations), determined by a characteristic material behavior specified by constitutive equations, which relate forces and deformations. Assuming only small displacements of the body, forces and deformations can be described by the Cauchy stress tensor σ and the infinitesimal strain tensor ε as defined in any textbook on continuum mechanics [47]. The characteristic property of a viscoelastic constitutive equation is that stress at a certain time instant depends on the entire history of strain and vice versa. Besides this memory hypothesis, three further assumptions lead to a quite general viscoelastic constitutive equation. These assumptions are • causality (no future stress resp. strain state can affect the current strain resp. stress state), • linearity (the principle of superposition holds) and • non-aging (the material behavior is independent of shifts in time).
The resulting viscoelastic constitutive equation is given by see [48][49][50]. Here, the functions J ijkl are known as creep functions and describe the strain increase of a material under an imposed unit step stress (given by the Heaviside function Θ(t)). Particularly, from (37) follows by integrating the Dirac delta function An alternative representation of the constitutive law in terms of stress relaxation can be obtained by reversing the roles of stress and strain in (37), see e.g., [48] The functions G ijkl are denoted as relaxation functions and represent the stress response to an imposed unit step strain, i.e., The general tensorial constitutive equations (37) or (39) can be simplified by using intrinsic symmetries of the stress and strain tensors as well as material symmetries. For the applications in this article, we consider the most simple case of isotropic material behavior. This leads to a description by only two remaining creep functions J h and J d such that where are the hydrostatic (subscript h) and deviatoric (subscript d) components of the strain tensor ε and the stress tensor σ, see [51] for a detailed derivation. Hence, (41) and (42) are functional equations with two scalar and independent creep functions J h and J d , that fully represent the isotropic constitutive law. An equivalent representation as (41) and (42) can be given in terms of two independent relaxation functions G h and G d . In the following, we will omit the subscripts h and d of the creep and relaxation functions when their general properties are discussed. From energy considerations and using the principle of fading memory [48] (Chapter 3), one can obtain monotonicity conditions for creep and relaxation functions and their first and second derivatives, i.e., Moreover, as we are interested in describing viscoelastic solids, we assume an initial elastic material response, a finite asymptotic creep and a non-zero asymptotic relaxation such that The most simple creep and relaxation functions that fulfill (44) and (45) are given by [48] (Chapter 1.5) where τ ε and τ σ are characteristic time quantities called retardation time and relaxation time, respectively. They determine the time scale, in which creep and relaxation take place, see Figure 2. As mechanistic idea, the creep function (46) and the relaxation function (47) are associated to a mechanical network model consisting of two springs (moduli E 0 , E 1 ) and a dashpot (viscosity η), named Zener model [52], see Figure 3. Connecting Hooke's law for a spring σ = Eε and Newton's law for a dashpot σ = ηε, the Zener model is represented by an ordinary differential equation [42]  It leads under the assumption of a unit step stress at t = 0 and the initial condition J(0) = 1 E 0 +E 1 (representing the instantaneous elasticity of the springs), to the creep function representation Moreover, assuming an imposed unit step strain at t = 0 and the initial elasticity G(0) = E 0 + E 1 , the relaxation function can be obtained from (48). The representations (49) and (50) are equivalent to (46) and (47), respectively, when defining the respective parameters by a comparison of coefficients.
The applicability of the Zener model to fit real material data is limited, as creep and relaxation processes usually proceed over a large time range. A possible improvement is given by the Prony series approach, considering several retardation and relaxation times τ ε,i , τ σ,i , i = 1, . . . , n, leading to creep and relaxation functions However, creep and relaxation functions of the form (51) and (52) can only fit the experimental creep and relaxation data in the measured time ranges and predict a fast decay of creep or a fast relaxation outside the measured time span, which results in a cumbersome extrapolation [1,2]. The most general representation of a viscoelastic constitutive law is obtained by assuming a continuum of retardation and relaxation times leading to retardation and relaxation spectra. For brevity, only the relaxation representation will be considered subsequently. The continuous version of (52) is given by where R σ denotes the relaxation spectrum in time, see [42,53]. An equivalent representation in terms of the frequency λ = 1 τ is given by where denotes the relaxation spectrum in frequency. A continuous relaxation spectrum is a characteristic property of fractional viscoelastic models and related to (12) as will be shown in the following section. Whereas creep and relaxation functions are quantities that represent viscoelastic material behavior under static deformation or loading conditions in the time domain, there are also quantities that describe the dynamic properties in the frequency domain, named the complex compliance and the complex modulus. These quantities determine the response of a linear viscoelastic model to an oscillatory input and will later be used to link the reformulated infinite state representation to fractional viscoelasticity. Particularly, consider the complex modulus E * , which relates the stressσ(t) to an oscillatory strain ε(t) = e iωt . Together with the one-dimensional version of (39) and the decomposition one obtains the representation in terms of the Laplace transform L ofĜ. Thereby, E * can be separated into real and imaginary part The quantities E * 1 and E * 2 are sometimes referred to as storage and loss modulus, respectively, see [48] (Chapter 1.6) and [50] (Chapter 3). This terminology becomes clear by using (58) in (57) The complex modulus can alternatively be given in terms of the relaxation spectrum. Substituting (54) in (58), interchanging integrals and using the Laplace transforms of sine and cosine, see [54], leads to as given in [53] (Chapter V.3).

Fractional Zener Model
Fractional calculus has been used by several authors in the thirties and forties of the last century as an empirical method to describe viscoelastic material behavior, see the references in [42] (Chapter 3). Particularly in [8], a fractional derivative was introduced to model intermediate material behavior interpolating between Hooke and Newton's law. Later in [13], fractional derivatives have been related to a new rheological element to be used in mechanical network models besides spring and dashpot, named springpot and depicted with a rhombus as symbol, see Figure 4. (Originally, the spelling "spring-pot" was used in [13]. We prefer an unhyphenated notation.) The springpot can be expressed in terms of different viscoelastic quantities. First, consider its slowly increasing creep function with parameters p > 0 and α ∈ (0, 1). Comparing (37) and (3), one can see as in [42] (Chapter 3), that (62) leads to a (one-dimensional) constitutive law which is in view of (5) equivalent to Moreover, the representations (39) and (64) reveal the associated relaxation function Using the infinite state representation (17) of (64), the relaxation spectrum of the springpot results with the help of (54) as S σ (λ) = pµ 1−α (λ).
The complex modulus of the springpot is obtained by substituting (66) in (60) and (61) such that the storage and loss moduli are given by where the expressions on the right-hand side result from Proposition 1. Hence, the storage and loss moduli of the springpot reappear in the reformulated infinite state representation (22). Particularly, the constitutive law (64) can be described by where E * 1 and E * 2 are assumed to fulfill (67) and (68). A springpot is the basic rheological element in fractional viscoelasticity. However, it is not sufficient to describe a viscoelastic solid, because (62) fulfills (44) but not (45). As generalization of the classical Zener model, one can consider a spring (modulus E 0 ) in parallel to a springpot (parameters p, α) with another spring (modulus E 1 ) in series, which is known as fractional Zener model [42], see Figure 4. The system is similar to (48) described by the functional equation A unit step stress at t = 0 and the initial condition J(0) = 1 , t ≥ 0, of the form (6). Applying the solution (10) of (6) to (71) leads to the creep function The relaxation function results from (70) similar to (47) as The creep behavior according to (72) is shown in Figure 5. The limit case α → 1 results in the classical Zener model (49), which shows creep on a small time range, whereas smaller values of α ∈ (0, 1) lead to an extended time window of creep. Comparing (46) with (72), one obtains a generalized retardation time which again determines the characteristic time scale on which creep occurs. Similarly, by comparison with (47), the relaxation time is given by More general fractional viscoelastic models can similar as (51) and (52) be formulated as

Finite Element Method
In order to describe the viscoelastic material behavior of an arbitrary structure by a fractional constitutive law, the above models have to be implemented in the finite element method. The formulation below, being the main result of this paper, uses the reformulated infinite state representation and the RISS scheme. It is generally applicable to (possibly nonlinear) fractional viscoelastic models. In the following, the method is demonstrated for a fractional Zener model.

Formulation of the Fractional Zener Model for a 3D Continuum
As starting point, consider once again a (one-dimensional) fractional Zener model ( Figure 6) as introduced in Section 3.2. An equivalent description of (70) is given by the equations where the internal variable ε f denotes the strain of the springpot in Figure 6. This representation complies with the state variables approach for classical viscoelastic constitutive laws introduced in [49] (Chapter VII) and [55] (Chapters 3 and 9) or the internal variable model in [35] (Chapter 10). In the formulation (78a) and (78b), the position dependence of stress and strain is explicitly taken into consideration. Thus, the model parameters are assumed to be independent of the position variable x, i.e., a homogeneous material is considered. Furthermore, substitution of the reformulated infinite state representation (22) of the fractional derivative of ε f in (78b) yields For the three-dimensional isotropic case, we assume a fractional Zener model (79a)-(79d) for both the hydrostatic and deviatoric components. To formulate the associated relations, we use the Cauchy stress tensor σ = σ xx σ yy σ zz σ xy σ yz σ zx T (80) and the linear strain tensorε = ε xx ε yy ε zz γ xy γ yz γ zx T (81) in Voigt notation. The separation in hydrostatic and deviatoric parts is, equivalent to (43), given byσ =σ h +σ d = T hσ + T dσ , with the matrices as introduced in [2]. Finally, adding up the hydrostatic and deviatoric versions of (79a)-(79d) results using (82) in the constitutive law

FEM Formulation
The equilibrium equations for a deformable body represented by a spatial domain Ω ⊂ R 3 can be formulated in terms of the displacement field being kinematically related to the linear strain ε. In case of a Cartesian coordinate system, the strain-displacement relation is given by The principle of virtual work including the virtual work of inertia δW dyn , of internal forces δW int and of external forces δW ext yields the equation which is valid for all smooth virtual displacement fields δu. Herein, ρ is the (volumetric) mass density in Ω, b is the volume density of external forces in Ω and t is the surface density of external forces on the boundary ∂Ω. For the displacement field u, we consider a finite element discretization where the matrix H contains the chosen shape functions and v the associated nodal displacements. Furthermore, the linear strain field is discretized as where the matrix B contains spatial derivatives of the shape functions used in H according to a strain-displacement relation such as (86). Furthermore, similar relations are assumed for the internal variablesε f ,Z andz in the form with generalized internal nodal displacement variables v f , Z and z. Using (89) and (90) in the principle of virtual work (88) results in Hence, inserting the constitutive law (84a) in (92) leads to the equations of motion where and the generalized system matrices have been introduced. Thereby, (94) and (95) are the classical mass and stiffness matrices, respectively and (97) is the generalized external force vector. In addition to these quantities, which are typical in FEM, the expression (96) represents a generalized system matrix expressing the influence of the internal states. Furthermore, multiplication of the differential equations (84b), (84c) and (84d) by B T , using (91) and integration leads to The equations of motion (93) together with (99a)-(99c) represent a reformulated system of FODEs, which can be solved numerically, given initial and boundary conditions. The states of the system include the nodal displacements v as well as internal states v f and infinite states Z, z.

Numerical Implementation
As (93) and (99a)-(99c) result from a fractional-order differential equation, the scheme RISS can be used to solve associated initial and and boundary value problems. Therefore, the Equations (99a)-(99c) are approximated using (33) as The integration of (93) and (100a)-(100c) can be performed as described in Section 2.4. The following fundamental benchmark examples demonstrate the application of the method and have closed form exact solutions, allowing for a proper error analysis.

Example 1.
Consider the one-dimensional example of a two-node rod element of length l and cross sectional area A that is fixed at one end and loaded by a constant force F at initial time t = 0 at the other end, see Figure 7. For this simple case, the quantities of stress (80) and strain (81) become scalar and the finite element scheme can be derived directly from (79a)-(79d). For the FEM discretization, we consider linear shape functions such that where v 1 and v 2 represent the displacement of left and right node, respectively. Neglecting inertia effects, (93) and (99a)-(99c) can be formulated for the given case as The numerical solution of (102a)-(102d) is carried out using the scheme RISS with quadrature parameters J = 10, K = 25, η 1 = 10 −5 , η K = 10 5 , As boundary conditions, the clamped left-hand side is considered, i.e., v 1 = 0, and the force F acting on the right-hand side of the rod. Furthermore, the rod is assumed to be fully relaxed initially such that all system states fulfill zero initial conditions.
The numerical solution can be evaluated by the closed form solution resulting from the creep function of the fractional Zener model. In view of the spatially constant unit-step stress field and the creep function J from (72), one obtains The relative error ∆ over time t between numerical and closed form solution for various values of α is shown in Figure 8. Therein, the Mittag-Leffler function has been computed as proposed in [56][57][58].

Example 2.
As a second benchmark problem, a quadratic plate is modeled as a two-dimensional finite element under plane strain conditions ( Figure 9). Accordingly, a reduced strain statẽ is considered resulting in a reduced stress statẽ Thereby, the stress σ zz in normal direction does not vanish but may be expressed depending on the plane stress variables in (105). A finite element discretization can be given by with displacements v 1x , . . . , v 4y in axial directions of the reference frame for the four nodes and bilinear shape functions Furthermore, in order to obtain a closed form solution of the problem, we consider boundary conditions such that a purely deviatoric deformation occurs. In particular, nodes 3 and 4 on the bottom are clamped, nodes 1 and 2 are fixed in y-direction and there is a constant distributed shear loading τ in x-direction on the face between nodes 1 and 2 applied as a step function at time t = 0. This leads to an external force vector with the components according to (97), Again, neglecting inertia, (93) and (99a)-(99c) can be formulated as where Q d is computed as in (98) with a reduced matrix Note that the parameters of the hydrostatic part of the constitutive law may be neglected due to the purely deviatoric deformation. Inserting the displacement boundary conditions v 1y = v 2y = v 3x = v 3y = v 4x = v 4y = 0, one can solve (110a)-(110d) for v 1x and v 2x . Moreover, a closed form solution can be obtained from the spatially constant stress boundary condition σ xy (x, 1, t) = τΘ(t), x ∈ (−1, 1) and a deviatoric creep function J d similar as in (72). The resulting nodal displacement is given by v 1x (t) = v 2x (t) = 2γ xy (x, 1, t) = 2τ J d (t) The relative error ∆ over time t between the numerical and the closed form solution is depicted in Figure 10.

Discussion
In Section 4, we have proposed and tested an implementation of the fractional Zener model in the finite element method using the reformulated infinite state scheme (RISS). The theoretical basis and motivation for fractional constitutive models is given in Section 3. In particular, it has been shown how the (reformulated) infinite state representation can be interpreted in terms of viscoelastic quantities of the springpot.
In contrast to most of the existing schemes, the method in Section 4 is not based on a Grünwald-Letnikov formulation but uses the infinite state scheme RISS. Accordingly, the histories of stress and strain magnitudes do not have to be stored and used in every computational step but are implicitly given in terms of the internal infinite states. The scheme, although given for a fractional Zener model, can be generalized for any other fractional constitutive law. The fractional Zener model has been chosen to evaluate the method, as it is the most simple fractional model of a viscoelastic solid with known creep and relaxation functions in closed form. So far, the method has been applied to simple benchmark problems for one-and two-dimensional finite elements, showing very small errors. An application to more complex problems and different element types and loading conditions are necessary in the future to further evaluate the scheme. Moreover, a comparison of the method to other common FEM formulations using infinite state schemes such as [38,39] should be performed.

Conflicts of Interest:
The authors declare no conflict of interest.