Higher order temporal finite element methods through mixed formalisms

The extended framework of Hamilton’s principle and the mixed convolved action principle provide new rigorous weak variational formalism for a broad range of initial boundary value problems in mathematical physics and mechanics. In this paper, their potential when adopting temporally higher order approximations is investigated. The classical single-degree-of-freedom dynamical systems are primarily considered to validate and to investigate the performance of the numerical algorithms developed from both formulations. For the undamped system, all the algorithms are symplectic and unconditionally stable with respect to the time step. For the damped system, they are shown to be accurate with good convergence characteristics.


Introduction
Despite of its origin in particle dynamics, Hamilton's principle (Hamilton 1834;Hamilton 1835) has been with us for a long time throughout broad range of mathematical physics (Bretherton 1970;Gossick 1967;Landau & Lifshitz 1975;Slawinski 2003;Tiersten 1967). However, it suffers from two main difficulties such as (i) use of end-point constraints and (ii) adoption of Rayleigh's dissipation for non-conservative systems. The first difficulty relates to the proper use of initial conditions resulting from the restrictions on the function variations. In Hamilton's principle, the variations vanish at the end points of the time interval, which, in turn, implies that the functions are known at these two instants. For a typical dynamic problem, one does not know how the considered system evolves at the end of the time interval. Usually, this is the main objective of the analysis, which means that there may be a serious philosophical or mathematical inconsistency in Hamilton's principle. Second difficulty relates to the inability to incorporate irreversible phenomena. Hamilton's principle itself only applies to conservative systems. With Rayleigh's dissipation (Rayleigh 1877), irreversible processes can be brought into the framework of Hamilton's principle. However, this approach is not satisfactory in a strict mathematical sense, since the variation of Rayleigh's dissipation enters in an ad-hoc manner.
Historically, to resolve such difficulties in Hamilton's principle, Tonti (Tonti 1973) suggested that convolution should replace the inner product for variational methods in initial value problems. Somewhat earlier, (Gurtin 1963;Gurtin 1964a,b) introduced the convolution functional, and could reduce the initial value problem to an equivalent boundary value problem. However, the functional by Gurtin is complicated and it never can recover the original strong form. Following the ideas of Tonti and Gurtin, Oden and Reddy (Oden and Reddy 1976) extended the formulation to a large class of initial boundary problems in mechanics, especially for Hellinger-Reissner type mixed principles.
Recently, two new variational frameworks for elastodynamics such as extended framework of Hamilton's principle (EHP, (Kim et al. 2013)) and mixed convolved action principle (MCAP, (Dargush and Kim 2012)) were established by using mixed variables. While EHP adopts a mixed Lagrangian formalism given in (Apostolakis and Dargush 2012;Apostolakis and Dargush 2013a;Sivaselvan et al. 2009;Sivaselvan and Reinhorn 2006), it provides a new and simple framework that correctly accounts for initial conditions within Hamilton's principle. EHP resides in an incomplete variational framework since it requires Rayleigh's function for dissipative systems and cannot define the functional action, explicitly. On the other hand, MCAP clearly resolves long-standing problems in Hamilton's principle. With MCAP, a single scalar functional action provides the governing differential equations, along with all the pertinent boundary and initial conditions for conservative and non-conservative linear systems. Thus, in theoretical aspects, MCAP is certainly preferred rather than EHP, however, there still remains a challenge for MCAP to have the generalized framework of other than linear problems. While EHP can be numerically implemented for viscoplasticity continuum dynamics, MCAP currently suffers from finding the explicit functional action for that problem. Since both methods provide sound basis to develop various space-time finite element methods for linear initial boundary value problems, here, the focus is initially on investigating their potential when employing higher-order temporal approximations.
The remainder of the paper is organized as follows. Next, in Section New variational formalisms, some relevant background on EHP and MCAP are provided, especially for the SDOF Kelvin-Voigt system. In Section Numerical implementation, discretization scheme and numerical algorithms are provided when temporally higher-order approximations are adopted in both approaches. Basic numerical properties of the developed methods are closely examined in Section Basic numerical properties. Then, some numerical examples are presented to investigate and to validate all of these developed algorithms for practical problems of the forced vibration in Section Numerical examples. Finally, some conclusions are provided in Section Conclusions.

New variational formalisms
In this section, new variational frameworks for the SDOF Kelvin-Voigt system displayed in Figure 1 were reviewed for the development of higher order temporal finite element methods from both approaches.
With mass m, damping coefficient c, the known applied forcef t ð Þ with time t, and stiffness k = 1/a with a representing the flexibility, EHP and MCAP could formulate the variational framework for this model in terms of the displacement of the mass u(t) and the impulse of the internal force J(t) in the spring.

Weak form for the Kelvin-Voigt model in EHP
Following the ideas in (Sivaselvan and Reinhorn 2006), the EHP associated with this problem defines Lagrangian L and Rayleigh's dissipation φ as where a superposed dot represents a derivative with respect to time. Then, the functional action A for the fixed time interval from t 0 to t is given by and, in EHP, the first variation of A is newly defined as by adding the counterparts (the underlined terms in Eq. (4)) to the terms without the end-point constraints in Hamilton's principle.
Such added terms have effect on confining a dynamical system to evolve uniquely from start to end with the unspecified values at the ends of the time interval such as_ u t 0 ð Þ;û t 0 ð Þ;_ u t ð Þ andû t ð Þ. Then, interpreting the unspecified initial terms as sequentially assigning the known initial values completes this formulation. Thus, in EHP, the given initial velocity _ u 0 is assigned first, and the given initial displacement u 0 is assigned next bŷ and The subsequent zero-valued term (6) needs not appear explicitly in the new action variation, so that the new definition (4) with the sequential assigning process such as (5) and (6) can properly account for the initial value problems. It should be noted that in EHP, the dependent initial condition J 0 can be identified by whereĵ 0 is the initial internal impulse of the known applied forcef given bŷ In Eq. (8), the time interval [−∞, t 0 ] is used to represent that this is the time interval before the initial time we are considering.
To check this, let us substitute Eqs. (1), (2) into Eq. (4). Then, we have Doing integration by parts on m _ u δ _ u; a _ J δ _ J ; and −u δ _ J in Eq. (9) yields For arbitrary variations of δu and δJ for the time interval [t 0 , t], the governing differential equations are given by along with constitutive relation as With the underlined terms in Eq. (10), the trajectory of the damped oscillator is firstly uniquely confined bŷ while the given initial conditions are identified sequentially by Eq. (5), Eq. (6) and Eq. (7). Thus, with EHP, Hamilton's principle can account for compatible initial conditions to the strong form. It is not a complete variational method, since it still requires the Rayleigh's dissipation for a non-conservative process and the first variation of the functional action cannot yield the proper weak form explicitly. However, the framework is quite simple and it can be readily applied to problems other than linear elasticity with the use of Rayleigh's dissipation.
For a representative example, let us consider SDOF elasto-viscoplastic model in Figure 2. Rayleigh's dissipation to define rate-deformation for the slider-dashpot _ u 1 can be given by in terms of Macaulay bracket 〈 ⋅ 〉 and absolute value of _ J whereby η and F y represent viscosity and yield force, respectively. Thus, in EHP, the action variation for this model is defined by adding up Eq. (4) and δA φ where the underlined term represents the rate-deformation for the slider-dashpot, _ u 1 . Note that the adding terms (16) are the counterparts to the terms without the end-point constraints in Hamilton's principle that obtained from the compatibility condition With Eq. (4) and Eqs. (14, 15 and 16), the governing differential equations for Figure 2.
are properly recovered in EHP along with proper initial conditions such as Eqs. (5, 6 and 7) andû 1 at t 0 .

Weak form for the Kelvin-Voigt model in MCAP
As well described in (Dargush and Kim 2012), MCAP defines the convolved action for the SDOF Kelvin-Voigt damped oscillator as where a superimposed arc represents a temporal left Riemann-Liouville semi-derivative. Referred to (Oldham and Spanier 1974;Samko et al. 1993), this is defined by where Γ(⋅) denotes the Gamma function. In Eq. (19), the symbol * represents the convolution of two functions over time, such that Meanwhile, the last termĵ 0 ð Þ in Eq. (19) represents the initial impulse corresponding tof t ð Þ that is given bŷ In MCAP, the stationarity of the action (19) yields the following weak form in time After performing classical and fractional integration by parts on the appropriate terms in Eq. (23) as follows (Apostolakis and Dargush 2012), we have For the sake of completeness, the fractional integration by parts formula is given For arbitrary variations of u and J, Eq. (24) emanates the governing differential equations in mixed forms as along with the proper initial conditions Note that the initial variations such as δu (0) and δJ(0) vanish due to Eq. (27). In other words, in MCAP, we can identify the dependent initial conditions such as J(0) andĵ 0 ð Þ from the usual given initial conditions u(0) and _ u 0 ð Þ as well as the known initial impulseĵ 0 ð Þ.
As shown in Eq. (24) and Eqs. (26, 27), every governing equations and initial conditions are satisfied weakly in MCAP, where it incorporates both conservative and non-conservative components within the unified functional action (19). Thus, it resolves the long-standing problem in Hamilton's principle. However, MCAP still requires a generalized framework to embrace various irreversible phenomena. In particular, currently, it does not have the functional action for the problem shown in Figure 2. Also, it should be noted that any pair of complementary order of fractional derivatives in Eq. (19) yields Eqs. (26, 27) due to the integration by parts property of complementary order of fractional derivatives for 0 < α < 1.

Numerical implementation
The weak form (9) in EHP and the weak form (23) in MCAP include, at most, first derivatives of the primary variables u(t) and J(t) as well as the variations δu(t) and δJ(t).
Consequently, we have C 0 temporal continuity requirement on primary variables and the variations, thus, there are many cases to develop higher order temporal finite element methods. As we shall see in this section, three kinds of quadratic temporal finite element methods in each framework are developed, since they are practically sufficient and accurate in computational aspects as discussed next. The numerical methods developed here are summarized in Table 1.

Algorithms from EHP
By introducing the fixed time step h for each time duration, that is, t r = r h, Eq. (9) can be written where δA r represents the action variation in the r th time duration [t r − 1 , t r ]. Also, p represents linear momentum, wherê For t r − 1 ≤ τ ≤ t r , temporally linear shape functions such as L r − 1 at t r − 1 and L r at t r are given by Also, by introducing the center point t c for the time interval [t r − 1 , t r ] as temporally quadratic shape functions Q r − 1 at t r − 1 , Q r at t r , and Q c at t c can be written as With linear temporal shape functions (31)-(32) and quadratic temporal shape functions (34)-(36), we can develop every algorithms of EHP presented in Table 1.
From either Eq. (49) or Eq. (50), we can express J c in terms of J r − 1 , J r , u r − 1 and u r . Then, replacing J c in the other independent equations with the equation of J c (J r − 1 , J r , u r − 1 , u r ) yields the matrix equation of Jquad algorithm as where X is given by Similarly, we have the Uquad algorithm as X h 2 2 6 6 6 6 6 4 3 7 7 7 7 7 5 u r p r J r 8 < : X h 2 2 6 6 6 6 6 4 3 7 7 7 7 7 5 where Y is given by Also, we have the UJquad algorithm as with the adequate substitution of u c and J c in terms of J r − 1 , J r , u r − 1 , and u r .

Algorithms from MCAP
Previously, MCAP was numerically implemented through linear temporal shape functions for classical SDOF oscillators and systems that utilize fractional-derivative constitutive models by (Dargush 2012). Here, continuing through this line, but, the quadratic temporal finite element methods are developed.
As well described in (Dargush 2012), for any non-negative integer m and n, we have the following relation for the convolution of the semi-derivatives of power functions.
To evaluate the convolution of semi-derivatives of polynomial shape functions, here, Eq. (56) is frequently used.

Since we cannot have summation form of the action variation in convolution integral
(that is, δA≠ X N r¼1 δA r ), let us consider the action variation over one time-step [0, h] as where temporally linear and quadratic shape functions of t (0 ≤ t ≤ h) are defined as Then, subsequent approximations are given by Now, let us consider Jquad algorithm for a representative one.
Then, by letting t → h in Eq. (69) due to the underlined term in Eq. (57), Eq. (69) yields In a similar way, and for the viscous dissipation term With evaluation of typical integer order convolution components in Eq. (57), we have the following discretized weak form of Jquad: or simply where Symplectic nature For the undamped case with no external forcing (conservative harmonic oscillator), Eqs. (87, 88 and 89) reduce to x n ¼A x n−1 ð91Þ where A left and A right in each algorithm are identified in Table 2 and Table 3. Notice that the magnitude of all the eigenvalues including complex conjugate pairs in each Table is exactly equal to 1, which can be written simply as Consequently, in addition to being time reversible, all the presented quadratic temporal finite element algorithms are also symplectic, energy conserving, and unconditionally stable for the undamped case.

Period elongation property in each method
To check the period elongation property in each developed method, the method by (Bathe 1996;Bathe and Wilson 1972) is used for free vibration of the undamped oscillator, where the ratio of the time-step h to the natural period T n is a control parameter. Also, Newmark's constant average acceleration method and Newmark's linear acceleration method are adopted for the references.
As shown in Figure 3, the numerical dispersion property from EHP and MCAP is exactly the same as Newmark's linear acceleration method, when either the primary variable u or J is quadratically approximated. On the other hand, when u and J are quadratically approximated, UJquad algorithm in each method has the same numerical dispersion property better than Newmark's linear acceleration method. Note that all while damping coefficient c = 0.2 π is fixed to deliver a non-dimensional damping ratio ξ = 0.05.
while the time step is fixed as h = 0.02. the developed methods are unconditionally stable, while Newmark's linear acceleration method is a conditionally stable algorithm with the criterion h/T n ≤ 0.551 . In computational aspects, compared to Newmark's constant average acceleration and Newmark's linear acceleration method, all the developed computational methods seem practically sufficient and accurate, since they have symplectic, unconditionally stable, and less or equivalent period elongation properties, and this is the main reason that only quadratic temporal finite element methods are developed here.

Numerical examples
For all of the numerical examples considered here, with no loss of generality, the model parameters are taken in non-dimensional form. In particular, let m = 1 and  a = 1/(4 π 2 ), thus, providing a natural period T n = 1 in the SDOF Kelvin-Voigt damped oscillator.
Two loading cases with zero initial conditions are considered for numerical simulation. The first one is an applied force in the formf t ð Þ ¼ f 0 sin ω 0 t ð Þ with f 0 = 100 and ω 0 = 10, and the other is 1940 El-Centro loading. The additional parameters for each loading case are summarized in Table 6.
For the references, the results obtained from each developed method are compared to an exact solution for the sinusoidal loading, while the results from Newmark's linear acceleration method in OpenSees (Mckeena et al. 2013;McKenna 2011) are additionally provided. For El-Centro loading, the results from each developed method are compared to those from Newmark's linear acceleration method in OpenSees.  As seen from the results, all the developed methods have better convergence characteristics compared to Newmark's linear acceleration methods under sinusoidal loading. In particular, UJquad algorithm in each framework shows the most accurate results.

Simulation results under 1940 El-Centro loading
The results from 1940 El-Centro loading analysis are displayed in Figures 11, 12, 13, 14, 15 and 16. In each figure, the Uquad and Jquad algorithms yield the exactly same results, while there are slight differences between the newly developed methods and  Newmark's linear acceleration method. In practical aspects, these differences seem negligible, but, note that all the developed methods are unconditionally stable that it may be advantageous to have the outlined results before the detailed analysis with the new methods.

Conclusions
In recent papers, through mixed formulation, two new variational frameworks such as EHP and MCAP were formulated for dynamical systems. Theoretically, MCAP is preferred to EHP, because unlike previous variational approaches, MCAP does not require any dissipation function with ad-hoc rules for taking variations, restrictions on the variations at the ends of the time interval, and external specification of initial conditions. However, there still remains a challenge for MCAP to have a generalized

Figure 11
Results from EHP algorithms for El-Centro loading analysis (1% damping ratio). framework embracing various irreversible phenomena. On the other hand, EHP has a relatively simple framework: the action variation is newly defined by adding the counterparts to the terms without the end-point constraints in Hamilton's principle, which confines a dynamical system to evolve uniquely from start to end. Interpreting these additional terms as sequentially assigning the known initial values completes this formulation. It should be noted that EHP is not a complete variational method, since it still requires the Rayleigh's dissipation for a non-conservative process and it cannot define the functional action explicitly. Since both mixed formalism provide a rigorous foundation to develop various temporal finite element methods for linear elasticity, in this paper, their potential when adopting temporally higher order approximations is investigated for the classical SDOF Kelvin-Voigt damped system. Figure 12 Results from MCAP algorithms for El-Centro loading analysis (1% damping ratio).

Figure 13
Results from EHP algorithms for El-Centro loading analysis (3% damping ratio).
With the consideration of computational aspects, three quadratic temporal finite element methods are essentially developed from each mixed formalism. All the developed methods are symplectic and unconditionally stable for the undamped conservative harmonic oscillator. Also, from period elongation property studies, it is checked that all the developed methods are equivalent or superior to Newmark's linear acceleration method that is conditionally stable. For damped forced vibrations, all the developed methods are shown to be robust and to be accurate with good convergence characteristics. It should be noted that since the new methods utilize mixed formulations, there exists an inherent disadvantage in a significant increase of the degrees of freedom against Newmark's methods when dealing with other than SDOF systems. However, this may be somewhat compensated by the general characteristics of a mixed Figure 14 Results from MCAP algorithms for El-Centro loading analysis (3% damping ratio).

Figure 15
Results from EHP algorithms for El-Centro loading analysis (5% damping ratio).
As the original Hamilton's principle has been adopted in various applications, the applicability of EHP and MCAP are quite broad, spanning many fields of mathematical physics and engineering. Future work will be directed toward development of a generalized framework of MCAP, and applications of both formalisms to various engineering problems, following the ideas in (Fried 1969;Hulbert 1992;Hulbert and Hughes 1990;Li and Wiberg 1996;Pitarresi and Manolis 1991;Bar-Yoseph 1989;Apostolakis and Dargush 2013b).