Mixed finite element solution of time-dependent problems

Available online at www.sciencedirect.com

Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678 www.elsevier.com/locate/cma

Mixed finite element solution of time-dependent problems J.A. Teixeira de Freitas Departamento de Engenharia Civil e Arquitectura, Instituto Superior Te´cnico, Technical University of Lisbon, Avenue Rovisco Pais, 1049-001 Lisboa, Portugal Received 23 November 2006; received in revised form 25 January 2008; accepted 18 February 2008 Available online 26 February 2008

Abstract A mixed formulation of the finite element method is used to establish a higher-order incremental method for the solution of secondorder/hyperbolic problems. The displacement, the velocity and, optionally, the acceleration fields are approximated independently in time using hierarchical bases. The time approximation criterion preserves hyperbolicity in the sense that it replaces the solution of hyperbolic problems by the solution of uncoupled Helmholtz-type elliptic problems, which can be subsequently solved using the alternative methods currently in use for discretization of the space dimension. The development of the time integration procedure, the characterization of its performance in terms of stability, accuracy and convergence are illustrated using a polynomial time basis. In order to stress the fact that the procedure can be implemented using alternative time bases, a wavelet system is used in the solution of nonlinear, parabolic and hyperbolic problems. The method is well-suited to parallel processing and to large time stepping. The extension of its application to the solution of other than linear second-order/hyperbolic problems is discussed. Ó 2008 Elsevier B.V. All rights reserved. Keywords: Time integration; Finite elements; Hyperbolic and higher-order problems

1. Introduction The time integration procedure presented in Ref. [3] is extended to the solution of higher-order problems. This procedure is designed to meet the objective of replacing non-periodic hyperbolic problems by equivalent Helmholtz-type elliptic problems. This objective was set to maintain in the time domain the relatively higher performance levels found for the Trefftz variant of the finite element method in the solution of quasi-static and of periodic dynamic problems solved in the frequency domain. The concept proposed by Trefftz [15], presented in 1926 as an alternative to the method suggested by Ritz [13], generalizes the classical method of solution of initial/boundary value problems. The approximation basis is still extracted from the set of solutions of the homogeneous problem but its weights are now determined by enforcing the boundary conditions in weak form. When this concept is used in

E-mail address: [email protected] 0045-7825/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2008.02.014

the framework of the finite element method, as pioneered by Jirousek [7], the approximation basis embodies the physics of the problem, which ensures higher levels of modelling performance, see Ref. [16]. They are well illustrated in the literature in the context of elliptic and time-dependent problems, the latter solved in the frequency domain, as the Fourier mapping preserves hyperbolicity, in the sense that it replaces the hyperbolic problem by equivalent Helmholtz-type elliptic problems. However, when the usual (typically implicit) time integration procedures are applied prior to the discretization in time, the fundamental properties of the problem are destroyed, yielding Trefftz approximation bases that weaken substantially the trademark features of this variant of the finite element method. Thus the motivation to develop the procedure described here, which is a direct extension of the time domain the technique described in Ref. [5] for the development of consistent hybrid and mixed finite element formulations and models. This procedure is, in essence, a modal decomposition technique implemented in the time domain, as opposed

3658

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

to the traditional method of modal decomposition in the space domain. These alternative approaches and their distinguishing features, in the context of the Trefftz method, are reported in Ref. [4], where they are related, also, with the Fourier method and the virtual impulse method of Tamma et al., e.g. Ref. [14]. This technique releases the limiting assumption of periodicity (or periodic extension) that hinders Fourier analysis and can be implemented on alternative time bases, as illustrated in Ref. [3] for polynomial, trigonometric, radial and wavelet systems. 2. Layout of the presentation The first part of the paper addresses basic aspects of the time integration procedure, as applied to second-order problems. A typical finite element approach is used in Section 3 to illustrate the selection and the implementation of the approximation criteria. It is applied to a scalar problem to focus the presentation on essential aspects. The procedure used to set up a high-order, naturally hierarchical time approximation basis yielding to fully uncoupled velocity and integration rules is summarized in Section 4. The extension of the time integration method to the solution of second-order systems is presented in Section 5. A general assessment of the time integration procedure is presented in the second part of the paper. The main characteristics of the method are recalled first in Section 6 to support a general assessment of its numerical implementation. Particular emphasis is placed on its suitability to parallel processing. The other issues that are addressed are the application of the time integration procedure to other than linear second-order problems, namely nonlinear, higherorder and hyperbolic problems. The coupling of the time integration procedure with modal decomposition in space is also commented in Section 6. This issue is addressed simply for the sake of completeness, as the time integration procedure used here does not rely on modal decomposition. The third part of the paper addresses the assessment of stability and accuracy. The criteria established in the literature are used in Section 7 to show that the integration scheme implemented on a complete, linearly independent polynomial basis of degree p is unconditionally stable and yields approximations of order O(X2p+2). The performance is illustrated with benchmark tests. The fourth and last part of the paper addresses the semidiscretization of parabolic/hyperbolic problems. Thus, the material presented in Section 8 is not essential to characterize the method in concept and in performance. Its relevance derives directly from the motivation to develop a non-periodic but still spectral decomposition procedure capable of preserving features essential to the implementation of time domain problems using the Trefftz concept. The advantages thus gained are illustrated with two numerical applications. The relative assessment of stability and accuracy results is particularly difficult in the present context, where a wide variety of time integration methods, and variants of the same basic method, exist. The option followed here is to

select the original versions of the Newmark and virtual impulse methods. Because the strengths and weaknesses of the trend-setting method proposed by Newmark [10] are so well known, its use here serves the single purpose of establishing a common background using a well-established standard. This is why the a-variants of the method are not used instead, and explains, also, why there is no intention of using the results obtained to judge upon the relative global performance. Comparison of low- and high-order integration methods is meaningful only in the context of a specific application. The virtual impulse method, overviewed and extensively assessed in Ref. [14], represents a distinct and, in many ways, embracing approach to time integration. Thus the option to include it also in the same background. The original version of the method (OVIP) is used instead of its generalized version (GVIP) simply because this latter version recovers the exact solution of the linear scalar problem that supports the stability and accuracy analyses. As for the Newmark standard, the issue is not to use the results reported here to extrapolate in terms of global performance. Similar reasons justify the use of trapezoidal rules to establish the background of the results reported in Section 8 on the solution of parabolic and hyperbolic test problems modelling the propagation of shock waves. 3. Outline of the time integration method Weighted residual approaches to time integration can be traced back to Zienkiewicz et al. [19] and Wood [18]. The basic aspects of this approach are introduced next using the linear, second-order scalar problem MaðtÞ þ CvðtÞ þ KuðtÞ ¼ F ðtÞ; uð0Þ ¼ u0 ;

ð1Þ ð2Þ

vð0Þ ¼ v0 ;

ð3Þ

where u(t) is the displacement, v(t) the velocity and a(t) the acceleration; M, C and K define the mass, damping and stiffness coefficients of the system, and F is the forcing load. 3.1. Primary approximation in time Problems (1)–(3) can be solved in the framework of the finite element method using alternative formulations that range from single- to multi-field approximation criteria. Typically, the period of analysis is divided into increments, 0 6 t 6 Dt, and the displacement field is approximated in each time interval, uðtÞ ¼

N X

T n ðtÞun ¼ TðtÞu;

0 6 t 6 Dt

ð4Þ

n¼1

where row-vector T(t), with dimension N, contains the time approximation functions, Tn(t), and vector u lists the corresponding weights, un. These generalized displacements are determined enforcing the governing Eq. (1) in a weak, Galerkin form

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

Z

Dt

T ðMa þ Cv þ Ku  F Þdt ¼ 0;

ð5Þ

0

where vector T* denotes the transpose conjugate of rowb T , as it may be convenient to define vector T(t), T ¼ T the time approximation basis in the complex space. Different strategies can be followed in the implementation of Eq. (5) under approximation (4). In essence, they depend on the techniques that are used to enforce the initial conditions (2) and (3) and the definitions of velocity and, eventually, acceleration: _ vðtÞ ¼ uðtÞ; aðtÞ ¼ v_ ðtÞ;

0 6 t 6 Dt;

ð6Þ

0 6 t 6 Dt:

ð7Þ

The stronger the requirements to satisfy locally these conditions, the stronger the constraints set a priori on the approximation basis. The option followed here is to release the time basis from all but the essential constraints on completeness and linear independence. The expected gains are in flexibility, in what concerns both the selection of approximation bases and the applicability of the time integration procedure to different classes of problems. 3.2. Dependent approximations in time Still in the context of the finite element method and in what concerns the local enforcement of the domain and initial conditions of the problem, the cost paid by releasing the constraints placed on the primary approximation (4) is to move from single- to multi-field (hybrid and/or mixed) formulations. A mixed formulation in the present context corresponds to approximate independently the velocity field and, eventually (see Section 5.4), the acceleration field: vðtÞ ¼ TðtÞv; aðtÞ ¼ TðtÞa;

0 6 t 6 Dt; 0 6 t 6 Dt:

ð8Þ ð9Þ

As for the displacement approximation, the weights vn and an listed in vectors v and a define generalized velocities and accelerations, with no immediate physical meaning. Consequent upon the option of limiting the constraints on the approximation basis to those that ensure basic convergence requirements, approximations (4), (8) and (9) may satisfy neither the initial conditions (2) and (3) of the problem, nor the initial condition on the acceleration, as determined by the equilibrium condition (1). In addition, approximations (4), (8) and (9) are not constrained, either, to satisfy locally the velocity and acceleration definitions (6) and (7). Instead, and very much in the way the mixed formulations of the finite element method are developed, these definitions are enforced in weak form, still in the sense of Galerkin Z Dt _ T ðv  uÞdt ¼ 0; ð10Þ

3659

to produce the velocity and acceleration integration rules, which express the velocity and acceleration weights, vn and an, in terms of the displacement weights, un. It is in this sense that approximations (8) and (9) are termed dependent on the primary approximation (4). 3.3. Algebraic solving system The implementation of approximations (4), (8) and (9) in the weak form (5) of the equilibrium condition (1) yields the following uncoupled system of algebraic equations: Mam þ Cvm þ Kum ¼ F m ; 1 6 m 6 N Z Dt ð12Þ DtF ¼ H1 T F dt 0

where H1 is the inverse of the following (positive-definite, Hermitian) matrix: Z Dt DtH ¼ T Tdt: ð13Þ 0

The computational cost of inverting (or factorizing) matrix H is marginal, as its dimension is the dimension of the time approximation basis, N (e.g., N = p + 1 for a polynomial approximation of degree p in the time interval Dt). The weak forms (10) and (11) of the initial conditions (2) and (3) are treated similarly. However, they are first integrated by parts in order to induce in the approximation the effect of the initial conditions (2) and (3): Z Dt Z Dt    T vdt ¼ T ðDtÞuðDtÞ  T ð0Þu0  ð14Þ T_  udt; 0

Z

Dt

T adt ¼ T ðDtÞvðDtÞ  T ð0Þv0 

0

Z

0 Dt

T_  vdt:

ð15Þ

0

The equivalent to this operation is present, also, either explicitly or implicitly, in the development of finite element formulations for elliptic problems, with the purpose of enforcing non-essential boundary conditions. Eqs. (14) and (15) are solved next for approximations (4), (8) and (9), to yield the following integration rules: Dtv ¼ H1 Gu  X0 u0 ; 1

Dta ¼ H Gv  X0 v0 ;

ð16Þ ð17Þ

where matrix G (in general, non-Hermitian) is defined by Z Dt G ¼ T ðDtÞTðDtÞ  ð18Þ T_  Tdt 0

and vector X0 weights the contribution of the initial conditions of the problem: X0 ¼ H1 T ð0Þ:

ð19Þ

The eventual computational usefulness of this procedure pivots on the possibility of uncoupling fully the velocity and acceleration integration rules (16) and (17), into form

0

Z

0

Dt

T ða  v_ Þdt ¼ 0

ð11Þ

Dtvn ¼ Xn un  X0n u0 ;

ð20Þ

Dtan ¼ Xn vn  X0n v0 :

ð21Þ

3660

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

This objective is attained whenever the following relation holds for a diagonal matrix X: G ¼ HX: ð22Þ Under this condition, system (12) remains uncoupled after eliminating the velocity and acceleration weights using the integration rules (20) and (21):  2  Xn Xn M þ C þ K un ¼ F 0n þ F n ; 1 6 n 6 N ; Dt2 Dt ð23Þ   X X0n 0n X0n 0 M þ C u0 þ Mv0 : Fm ¼ Dt Dt Dt System (23) is solved for the generalized displacements, un, which are used to recover the velocity and acceleration weights through Eqs. (20) and (21) and obtain the approximate solutions (4), (8) and (9) to the scalar second-order problem (1)–(3). It is recalled that, in general, this solution will not recover the initial conditions set for the time increment, that is x0 6¼ Tð0Þx;

ð24Þ

where x is the vector that collects the displacement, velocity or acceleration weights, and x0 is the initial condition set for the same variable, u0, v0 or a0, respectively. The information presented next, namely the construction of time approximation bases satisfying the uncoupling condition (22) and the extension of the time integration procedure to the solution of second-order systems of equations, is used in Section 6 to assess the application of this procedure to other than linear second-order problems.

Step 1: Choose the family of functions to be used in the approximation, u(s), and define the dimension of the basis, N. Step 2: Set up (the symmetric/Hermitian) matrix, Z 1 h¼ u u ds 0

and compute its (real) eigenvalues, K, and (orthonormal) eigenvectors, E: hE ¼ EK:

ð26Þ

Step 3: Set up (the asymmetric/non-Hermitian) matrix, where u_ is the s-derivative of u: Z 1 g ¼ u ð1Þuð1Þ  u_  u ds: 0

Step 4: Set up the auxiliary matrix 1

T

1

g ¼ ðEK2 Þ gðEK2 Þ

ð27Þ

and compute its (complex conjugate) eigenvalues, X, and eigenvectors, E*: g E ¼ E X:

ð28Þ

Step 5: Define the linear combination matrix as follows: 1

Z ¼ EK2 E : This result, together with conditions (26)–(28) and definitions, H ¼ Z hZ;

4. Construction of the time approximation basis

G ¼ Z gZ

The procedure used to set up time approximation bases is presented first for non-periodic approximations. The only purpose of its subsequent specialization to periodic approximations is to clarify its relation with well-established spectral analysis techniques.

obtained inserting definition (25) in Eqs. (13) and (18), can be used to show that the uncoupling condition (22) holds for the time approximation basis.

4.1. Definitions Eq. (18) shows that, in general, matrix G is non-Hermitian, meaning that the uncoupling condition (22) cannot be secured using a standard diagonalization method. For this reason, the time approximation basis is written in form TðtÞ ¼ uðsÞZ; ð25Þ where row-vector u(s), with s = t/Dt, collects the (eventually complex) approximation functions, and matrix Z defines the linear combination that ensures uncoupling. It is assumed throughout that basis u(s), defined on the unit interval, is complete and linearly independent.

It is noted that the eigenvalue problem (26) is trivial when the basis is orthonormal. Moreover, the storage required to register the basic information needed to implement the time basis is limited to the eigenvalues Xn, vector X0 and matrices H1 and Z, each with the dimension of the time basis, N. A complete polynomial (orthonormal Legendre) basis, with degree N  1, is used in all applications reported below, with the exception of the nonlinear, parabolic and hyperbolic testing problems, which are solved using the wavelet system defined on the unit interval proposed by Daubechies [1] and Monasse and Perrier [9], as implemented in Ref. [12]. 4.3. Periodic bases

4.2. Solution procedure The procedure used to determine matrix Z, and, in consequence, the constants present in the integration rules (20) and (21), can be summarized as follows:

If the solution is assumed to be periodic over the time interval Dt, u(Dt) = u(0) and v(Dt) = v(0), and a periodic basis is used, T(Dt) = T(0), the boundary terms in Eqs. (14) and (15) cancel out, definition (18) simplifies to form

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

G¼

Z

Dt

T_  T dt

0

and the integration rules (20) and (21) become independent of the initial conditions, as expected: Dtvn = Xnun and Dtan = Xnvn. The procedure described above can be used to set up the time basis (25), using now a family u(s) of periodic functions, u(1) = u(0). In particular, it can be readily verified that the eigenvalues Xn identify with the Fourier spectral frequencies, Xn = 2npj, where j is the imaginary unit, when a periodic trigonometric basis is used in Eq. (25), with un(s) = exp(2npjs) and Z = I, the identity matrix. It is in this sense that the solving system (23) can be interpreted as the description of problem (1) in the spectrum of (non-dimensional) frequencies Xn generated by the approximation basis, and corrected by a forcing term determined by the initial conditions of the non-periodic problem.

The dual transformations of approximations (32)–(34) are still used to enforce on average the governing system (29) and the extension of the velocity and acceleration definitions (6) and (7) to each degree-of-freedom of the system Z Dt Z Dt Tb n ðMa þ Cv þ KuÞdt ¼ Tb n Fdt; ð35Þ 0

Z

0 Dt

_ Tb n ðv  uÞdt ¼ 0;

ð36Þ

Tb n ða  v_ Þdt ¼ 0;

ð37Þ

0

Z

Dt

0

where Tb n is the complex conjugate of Tn and 1 6 n 6 N. It is noted that the following relations hold, in consequence of constraint (22) on the complete time basis Z Dt Xn T n dt ¼ DtT n ðDtÞ; 0 N X

5. Solution of linear second-order systems

MaðtÞ þ CvðtÞ þ KuðtÞ ¼ FðtÞ; uð0Þ ¼ u0 ;

ð29Þ ð30Þ

vð0Þ ¼ v0 ;

ð31Þ

where u, v and a are the displacement, velocity and acceleration vectors, respectively. The system mass, damping and stiffness matrices, M, C and K, are assumed to be constant but not necessarily symmetric or Hermitian. The definition of the forcing load vector, F, is constrained only to be supported by the (analytical or numerical) information necessary to implement the integral expressions defined below in the text. 5.1. Discretization criteria Similarly to the scalar problem (1)–(3), assume that each displacement, velocity and acceleration degree-of-freedom, um(t), vm(t) and am(t), with 1 6 m 6 Nsys, is approximated in the time interval 0 6 t 6 Dt. Approximations (4), (8) and (9) generalize into form uðtÞ ¼

N X

T n ðtÞun ;

ð32Þ

T n ðtÞvn ;

ð33Þ

T n ðtÞan ;

ð34Þ

n¼1

vðtÞ ¼

N X n¼1

aðtÞ ¼

N X n¼1

where Tn(t) defines the nth term of a time approximation basis determined according to the conditions and the procedure defined in Section 4. Its dimension is, in general, substantially smaller than the dimension of system (29), N  Nsys.

Z

X0n

Dt

T n dt ¼ Dt

0

n¼1

Consider now the solution of the second-order system of equations, with dimension Nsys

3661

as required by appropriate combinations of conditions (36) and (37), which imply a null error on the average velocity and acceleration over the increment Dt: Z Dt vdt  uðDtÞ þ u0 ¼ 0; 0

Z

Dt

adt  vðDtÞ þ v0 ¼ 0:

0

5.2. Algebraic solving system The algebraic solving system is obtained processing Eqs. (35)–(37) in the manner described in Section 3 for their scalar counterparts. This procedure is briefly recalled next. The velocity and acceleration consistency conditions (36) and (37) are integrated by parts to enforce the initial conditions (30) and (31) and approximations (32)–(34) are enforced in the resulting expressions, which are then integrated in time. The following velocity and acceleration integration rules are obtained after enforcing the uncoupling condition (22), Dtvn ¼ Xn un  X0n u0 ;

ð38Þ

Dtan ¼ Xn vn  X0n v0 ;

ð39Þ

where Xn is the solution of the eigenvalue of problem (28) and X0n is defined by Eq. (19). As for the weak form (5) of the scalar problem, approximations (32)–(34) are inserted in Eq. (35) and integration in time is implemented to produce a set of N systems Z Dt N X Tb m F dt; Dt H mn ðMan þ Cvn þ Kun Þ ¼ n¼1

0

where Hmn is the general coefficient of matrix H, defined by Eq. (13). The H1-combination of these systems, together

3662

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

with the integration rules (38) and (39), yields the following sequence of N uncoupled algebraic solving systems:  2  Xn Xn C þ K un ¼ F0n þ Fn ; 1 6 n 6 N : M þ ð40Þ Dt2 Dt Vectors F0n and Fn define the influence of the initial conditions and of the forcing load,   X0n Xn X0n M þ C u0 þ Mv0 ; ð41Þ F0n ¼ Dt Dt Dt Z Dt N X H 1 nm Tb m F dt; ð42Þ Fn ¼ Dt 0 m¼1 where H 1 nm denotes the general coefficient of the Hermitian matrix H1. 5.3. Velocity and acceleration estimates The displacement vector estimate at instant t is determined from the primary approximation (32), combining the N displacement vector weights determined from the solution of system (40). Similarly, the velocity and acceleration estimates are obtained either applying sequentially the integration rules (38) and (39) and approximations (33) and (34), or using directly the expressions that result from this operation: N X DtvðtÞ ¼ T n ðtÞXn un  h0 ðtÞu0 ; ð43Þ n¼1

Dt2 aðtÞ ¼

N X

T n ðtÞX2n un  h1 ðtÞu0  h0 ðtÞDtv0 ;

ð44Þ

Step 1: Implement the procedure described in Section 4 to set up and store the necessary data on the supporting time approximation basis, with dimension N; Step 2: Set the process time frame, t0 ¼ t00 þ t, where 0 6 t 6 Dt is the time frame adopted in each increment; Step 3: For each term of the time approximation basis, 1 6 n 6 N, compute vectors (41) and (42) and set up and solve system (40); Step 4: Update the initial conditions u0 = u(Dt) and v0 = v(Dt) using Eqs. (32), (33) and (38), or Eq. (43); Step 5: Update the time frame, t00 ¼ t00 þ Dt, and return to Step 2 if t00 < t0max . Step 4 can be extended to include the computation of the acceleration vector estimate using Eqs. (34) and (39), or the alternative expression (44). For large time stepping, it can also be extended to determine the solution estimates at particular instants of each time interval. The procedure summarized above does not imply the use of constant time steps, Dt, and its specialization to the solution of first-order systems is immediate. 6. General assessment The time integration procedure presented above is assessed in this section in terms of numerical implementation and of its application to the solution of other than linear second-order problems. In order to support this assessment, the basic characteristics of the time integration procedure are recalled first.

n¼1

hk ðtÞ ¼

N X

Tn ðtÞXkn X0n :

n¼1

6.1. Basic characteristics The following statements result directly from the information presented in Sections 4 and 5:

5.4. Alternative acceleration estimates The analysis of the procedure summarized above shows that approximation (34) may not be used, as the integration by parts of Eq. (35), together with the consistency condition (38) on the displacement and velocity approximation weights, yields the same algebraic solving system (40) and, consequently, the same displacement and velocity estimates. However, the two alternatives that are then left to obtain the acceleration estimate, whenever deemed relevant to characterize the response of the system, are not as competitive. They are the local enforcement of the acceleration definition (7), written for each acceleration component and using the velocity estimate (43), or the local enforcement of the equilibrium condition, stated by system (29), for the displacement and velocity estimates. 5.5. Solution procedure The numerical implementation of the time integration procedure presented above is straightforward, and can be summarized as follows:

(S1) Approximation basis: The implementation of the time integration procedure can be applied using any complete, linearly independent basis, u(s) in Eq. (25). In particular (a) The implementation of the uncoupling condition (22) involves the solution of a (generally, trivial) eigenvalue problem with the dimension of the basis, N, independently of the type or the dimension of the problem to be discretized in time; (b) The derivation of the data that supports the implementation of the time basis, namely the eigenvalues, Xn, the initial solution corrective term, X0n, the spectral decomposition matrices, Z and H1, has a marginal cost, both in computational and storage terms, and is independent of the (eventually variable) time increment, Dt; (c) The approximation basis, u(s), is independent of the problem to be solved, but can be chosen to model particular forms of response, namely periodic, or periodically extended problems (in which case a Fou-

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

rier basis yields the spectral decomposition for frequency domain analyses), and problems localized both in time and in frequency (in which case the use of systems of wavelets is particularly well-suited). (S2) Approximation in time: The resulting time approximation basis, T(t) in Eqs. (32)–(34), is independent of the problem to be solved in two essential aspects: (a) The velocity and acceleration integration rules (38) and (39) are valid for both linear and nonlinear problems; (b) The dimension of the supporting eigenvalue problem, N, is, in general, substantially smaller than the dimension, Nsys, of the time-dependent problem to be solved. (S3) Approximation criteria: The implementation of the time discretization criteria, defined by Eqs. (32)–(37) shows that: (a) No constraints are placed on the definition of the forcing load vector, F, and the system mass, damping and stiffness matrices, M, C and K, which may not be symmetric or Hermitian; (b) Modal decomposition of the homogeneous system (29), with dimension Nsys, is called upon at no stage of the derivation of its discretized form, defined by Eq. (40); (c) System (40) is fully uncoupled, in the sense that each term in the set 1 6 n 6 N can be solved independently, with N  Nsys, in general; (d) When complete time bases are used (N even), system (40) has pairs of complex conjugate solutions, meaning that it suffices to solve either N/2 complex systems or N real systems. Although statement (S3) is written for the second-order problem analyzed so far, the comments below show that it remains valid for both higher-order and hyperbolic problems. 6.2. Numerical implementation of linear second-order systems The computational cost of the implementation of the procedure presented in Section 5 centres fully on the solution of system (40). The cost of all other operations involved is marginal, namely the construction of the time approximation basis, the computation of the forcing load vector and the updating of the finite element estimates (which can be based on techniques equivalent to fast Fourier transforms). The selection of a particular solution strategy must attend to the features of the application in question, while weighing the characteristics of the time discretization procedure summarized above. The first option open is to exploit the modal decomposition of system (40). The strengths and weaknesses of this

3663

option are well-established. Thus, and recalling statement (S3b), the objective of the application summarized in Section 6.3 is two-fold: to contrast modal decomposition with the approach followed here, and to illustrate the combination of the two approaches. The second option is to factorize the system matrix, which may be justifiable for weak approximations in time (time bases with small dimension, N) inherently associated with strategies based on (relatively) small time increments. The third option, which is known to be particularly competitive when applied to well-conditioned, highly-sparse systems, is to use iterative solvers, adopting, eventually, the current estimate, un, as the solution initializer. Independently of the option followed in the selection of the solver, the advantages offered by the uncoupling condition, recalled in statement (S3c), are essential to the efficient solution of large systems in parallel processing. In spite of their limitations, the numerical tests reported below show that the time integration method is competitive for large time stepping, particularly when implemented in parallel processing, for which it is naturally well-suited. As a rule, its implementation under small time stepping is costlier than that of well-tuned implicit methods, even when weak approximation bases (N small) are used. 6.3. Application with modal decomposition Under the relatively stronger constraints on linearity, symmetry and proportional damping, see statement (S3a), the dynamic eigenvalue problem K/n ¼ x2n M/n has real natural frequencies, xn, and real, orthogonal vibration modes, /n, to yield: /Tm M/n ¼ M n dmn ; /Tm K/n ¼ K n dmn ; /Tm C/n ¼ C n dmn ; where K n ¼ x2n M n , Cn = 2nnxnMn and dmn is the Kronecker symbol. The orthogonality conditions and the approximation in time uðtÞ ¼

N sys X

/n un ðtÞ

n¼1

are used to uncouple problem (29)–(31) in a set of Nsys scalar problems, (if all modes are used in the approximation); €un þ 2nn xn u_ n þ x2n un ¼ M 1 n fn ; T un ð0Þ ¼ u0n ¼ M 1 n /n Muo ; T vn ð0Þ ¼ v0n ¼ M 1 n /n Mvo ;

ð45Þ

fn ðtÞ ¼ /Tn FðtÞ: System (45) has known homogeneous and particular solutions, which are given in Appendix 1 for later use. This feature is recalled at this moment to stress that the solution of

3664

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

boundary/initial-value problems using the homogeneous solution of the governing equation, that is, Eq. (45) with fn = 0 in the present context, is the root of classical integration methods and central to the generalization suggested by Trefftz [15]. The same concept supports the development of the GVIP method presented and assessed in detail in Ref. [14]. It is in this strict sense that the GVIP method can be classified as a Trefftz-like method. Modal decomposition is not used here, as the uncoupling condition is relaxed and limited to the spectral decomposition of the time approximation basis. However, it is possible to call upon it and implement approximations (32)–(34) with

C, and K, eventually computed at each time increment, or regularly updated over the analysis process, as more adequate. Results (40) and (41) still hold, and definition (42) for the forcing load contribution is extended to account for the effect of the source of nonlinearity: Z Dt N X H 1 nm Tb m ðF þ FL  FNL Þdt: Fn ¼ Dt 0 m¼1

un ¼ /n un ; vn ¼ /n vn ;

6.5. Application to higher-order problems

an ¼ /n an ; pffiffiffiffiffiffiffiffiffi T n ðtÞ ¼ an Dt expððnn xn þ jxdn ÞtÞ;

ð46Þ

where un, vn and an are now (measures of) displacement, velocity and acceleration modal amplitudes, xdn is the damped frequency of the nth mode (see Appendix 1), and parameter an is so chosen as to ensure that matrix (13) has unit diagonal coefficients, Hnn = 1. After enforcing the modal orthogonality conditions, the uncoupled velocity and acceleration integration rules (38) and (39) are recovered in form Dtvn ¼ Gnn un  u0n ; Dtan ¼ Gnn vn  v0n ; where Gnn is the (real) diagonal coefficient of the matrix defined by Eq. (18). The scalar form thus found for the solving system (40) is  2  Gnn Gnn C n þ K n un ¼ F 0n þ F n 1 6 n 6 N sys ; Mn þ Dt2 Dt   Gnn 0 M n þ C n u0n þ M n v0n ; DtF n ¼ Dt Z Dt Tb n fn dt: DtF n ¼ 0

Naturally, the solution procedure summarized in Section 5 has to be extended to include an iteration cycle, typically based on Newton–Raphson or arc-length techniques.

The extension of the time integration procedure described here to the solution of higher-order problems is straightforward. Let the second-order problem (29)–(31) be extended to order ‘, ‘ X

Sk uðkÞ ðtÞ ¼ FðtÞ;

k¼0 ðkÞ

uðkÞ ð0Þ ¼ u0 ;

k ¼ 0; 1; . . . ; ‘  1;

where u(0)(t) is the displacement vector, u(t), and u(k)(t) denotes its kth time derivative. Approximations (32)–(34) and integration rules (38) and (39) generalize to uðkÞ ðtÞ ¼

N X

T n ðtÞuðkÞ n ;

k ¼ 0; 1; . . . ; ‘

n¼1 ðk1Þ

ðk1Þ DtuðkÞ  X0n u0 n ¼ Xn un

k ¼ 1; 2; . . . ; ‘

(with X0n = 0 for periodic problems) and the solving system (40) extends as follows for the nth-order solution, where result (42) still holds k ! ‘  X Xn Sk unð0Þ ¼ F0n þ Fn ; Dt k¼0  k‘1 ‘ ‘2 X ‘ X0n X X0n X Xn ðk1Þ ðkÞ 0 Fn ¼ Sk u0 þ Skp u0 : Dt k¼1 Dt p¼0 k¼pþ2 Dt

As before, the solution of the system above is defined by pairs of complex conjugate values, as the time approximation basis (46) has to be extended (for completeness) to include the complex conjugate terms, Tb n ðtÞ.

Thus, the procedure summarized in Section 5 can be readily adapted to the solution of higher-order problems.

6.4. Application to nonlinear problems

Hyperbolic problems are usually discretized by separation of variables in space and in time, which raises the question on the sequence that should be followed in the implementation of the semi-discretizations in time and in space. Obviously, the sequence used in their implementation is irrelevant when the approximation bases in space and in time are chosen to be problem-independent. This is typically the case of most of the (single- or multi-field) applications that have been reported, based on the different

Consequent upon statement (S2a), the effects of sources of nonlinearity are restricted to the average enforcement of the governing system, as stated by Eq. (35) for secondorder problems. Thus, the implementation of the time integration procedure may follow the usual practice of collecting in vector FNL the nonlinear internal forces and in vector FL = Ma + Cv + K u the internal force estimate obtained with linearized mass, damping and stiffness matrices, M,

6.6. Application to hyperbolic problems

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

methods of space discretization used presently. This option leads to the second-order problem (29)–(31), which can be solved using the many alternative time integration procedures currently available. They may also be solved using the procedure summarized in Section 5, which justifies the brief comment presented in Section 8. However, and still in the context of separation of variables in time and in space, when one of the bases is chosen to be problem-dependent, the first (semi-) discretization should be implemented on the problem-independent basis. A typical illustration is the modal decomposition approach summarized above, which chooses a problem-dependent time approximation basis: discretization in space has to be implemented first to set up the supporting time-dependent system (29)–(31). If modal decomposition in the space dimension of the problem is to be avoided, it becomes necessary to discretize first the time dimension of the hyperbolic problem. The resulting elliptic (Helmholtz-type) problem is then discretized in space and its formal solutions are used to set up the space approximation basis (thus its problem-dependent nature). This is the option that, in the context of the Trefftz variant of the finite element method, has motivated the development of the time integration procedure reported here. Its implementation is presented and illustrated in Section 8, after assessing in the following section the stability and accuracy of the supporting time integration method.

3665

where q(A) is the spectral radius (the largest eigenvalue module of A): qðAÞ ¼ expðnXÞ:

ð52Þ

7.1. Measures of stability and accuracy To assess the stability and the accuracy of a given time integration method, the estimate of the solution found for the same initial condition and forcing load is expressed in form (48) xan ¼ Aa xn1 þ La ; which is used to define, by analogy, e.g. Pegon [11], the numerical spectral radius as the largest eigenvalue module of the approximate amplification matrix qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðAÞ ¼ a2 þ b ð53Þ q and the numerical estimates for the (normalized) damped and undamped frequencies and for the damping ratio:  aÞ;  d ¼ tan1 ðb= X qffiffiffiffiffiffiffiffiffiffiffiffiffi  ¼X  d = 1  n2 ; X n ¼  ln q  ðAÞ=X:

ð54Þ

ð55Þ

Unconditional stability is ensured by condition, ðAÞ 6 1 q

ð56Þ

7. Stability and accuracy The assessment of the accuracy and stability of secondorder time integration procedures is implemented on the scalar problem (1)–(3), written in the canonical form 2

Dt2 a þ 2nXDtv þ X2 u ¼

X F; K

ð47Þ

and accuracy can be measured at the end of the time step, tn = tn1 + Dt, by the truncation error found on the nonscaled solution, sT ¼ f u v a g: en ¼ sn  san

ð57Þ

This assessment is complemented with the following measures on numerical distortion and dissipation:

where X ¼ xDt is the normalized angular frequency, X2 = Dt2K/M, and n is the damping ratio, 2nX = Dt C/ M. The exact solution (see Appendix 1) is used to define the (scaled) solution xTn ¼ fun Dtvn Dt2 an g at instant tn = tn1 + Dt in terms of the effects of the initial solution, xn1 at instant tn1, and of the forcing load:

 d =Xd ; eX ¼ 1  X  en ¼ 1  n=n:

xn ¼ Axn1 þ L:

According to the results obtained in Section 3 for the scalar second-order problem, and letting

ð48Þ

ð58Þ ð59Þ

7.2. Approximate solution

The definitions thus found for the forcing load vector, L, and for the amplification matrix, A, are recalled in Appendix 2, where Xd defines the normalized damped frequency: qffiffiffiffiffiffiffiffiffiffiffiffiffi Xd ¼ 1  n2 X:

the general expression for the (non-scaled) approximate solution is,

The exact amplification matrix has one null eigenvalue and a pair of conjugate eigenvalues,

sT ðtÞ ¼ f uðtÞ

k1 ¼ 0; k2 ¼ a þ bj; a ¼ qðAÞ cosðXd Þ; b ¼ qðAÞ sinðXd Þ;

k3 ¼ a  bj;

ð49Þ ð50Þ ð51Þ

fm ðtÞ ¼ T m ðtÞDm ðF m þ F 0m Þ;

ð60Þ

2 2 D1 m ¼ 1 þ 2nX=Xm þ X =Xm

vðtÞ aðtÞ g ¼

N  X2 X fm ðtÞ X2 m K n¼1

X1 m

1



to yield the following definitions for the amplification matrix and for the forcing load vector:

3666

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

2

c00 þ 2nXc01 6 Aa ¼ 4 X2 c01

c01 c00

X2 c00

2nXc00  X2 c01

3 0 7 0 5; 0

LTa ¼ fLa2 La1 La0 g; N X2 X T n ðDtÞF n Dn =Xmn ; Lam ¼ K n¼1 ckm

¼

N X

7.4. Stability ð61Þ

ð62Þ

N X

ð63Þ

k H 1 mn an ;

ð64Þ

n¼1

akm ¼ ðk!Dtkþ1 Þ

1

Z

According to the results presented in Appendix 3, see Eq. (3.1), the definitions above can be written as follows for a polynomial basis with degree p, a ¼

Dt

T^ m ðtÞtk dt:

a ¼ c00 þ nXc01 ;  ¼ X d c0 : b 1

T n ðDtÞbkn Dn =Xmn ;

n¼1

bkm ¼

The approximate amplification matrix (61) presents the same eigenvalue structure (49) of the exact solution, to yield:

ð65Þ ¼ b

where the coefficients of matrix a and vector b are constants proportional to residues on the normalized frequency of order O(X2p+2) and O(X2p+2k), respectively. The emulation of multi-step methods has long been used to assess (and improve) the accuracy of single-step methods of integration, e.g. Ref. [6]. It can be verified that the multistep equation for the exact problem unþ1 ¼ A1 un  A2 un1 þ ðL2 Þnþ1  ðA22 L2  A12 L1 Þn ; where A1 = A11 + A22 is the trace and A2 = A11A22  A21A12 is the (non-trivial) principal minor of the amplification matrix, takes the same form for the approximate solution estimates defined by Eqs. (61) and (62). The orders of approximation defined above are recovered when the multi-step equation is used to estimate the truncation error of the solution at instant tn+1 = tn + Dt, enforcing the displacement estimate obtained at instant tn = tn1 + Dt and the exact solution defined at instant tn1.

Xd

d n Xn þ OðX2pþ2 Þ; ðn þ 1Þ!

d n Xn þ OðX2pþ2 Þ; ðn þ 1Þ!

ð66Þ ð67Þ

where parameters dn are the damping ratio dependent functions defined by Eq. (3.1) in the same appendix. These coefficients can be used to prove the convergence of the numerical solution to the exact eigenvalue solutions (50) and (51) expressed in power series form: a¼

1 X

ðn þ 1 þ nXÞ

n¼0

b¼

en ¼ asn1 þ bF ;

2pþ1 X n¼0

7.3. Truncation error The power series expansions of the exact and approximate amplification matrices are defined in Appendix 4, and the expansions obtained for the exact and approximate definitions of the forcing load are given in Appendix 5 for a polynomial load of degree k, and intensity F . These results can be used to show that the truncation error (57) on the non-scaled variables, namely on the displacement, velocity and acceleration estimates, can be defined in form,

ðn þ 1 þ nXÞ

n¼0

0

The results above show that the discriminant defined by Eq. (60) and the (double-index) parameters defined by Eqs. (63)–(65) play an important role in the derivation of the results presented here. Their properties and the information necessary to define their power series expansions on the normalized angular frequency, X, are summarized in Appendix 3 for the complete, linearly independent, polynomial basis, u(s), with degree p = N  1, see Eq. (25), that is used in all linear tests reported in this section.

2pþ1 X

1 X n¼0

Xd

d n Xn ; ðn þ 1Þ!

d n Xn : ðn þ 1Þ!

In addition, results (66) and (67) can be used, together with definition (52), to establish the following expression for the numerical spectral radius definition (53): ðAÞ ¼ qðAÞ  OðX2pþ2 Þ: q

ð68Þ

The values obtained for the numerical spectral radius implementing definition (53) with the same polynomial time basis with dimension N (degree N  1), are shown in Fig. 1 for different damping ratios. They are compared with the values obtained with the (trapezoidal) Newmark (NMK) scheme [10] and with the original version of the virtual impulse (OVIP) method [14]. The improved performance is evident in the high frequency range, even when weak time bases are used, namely with the linear approximation, N = 2. These results show, also, that the unconditional stability condition (56) holds for all ranges of frequency and (under) damping. The exact value for the spectral radius is not recovered with a null damping ratio because the (positive) residue in Eq. (68) does not tend to zero with the damping ratio: ðAÞ ¼ 1  OðX2pþ2 Þ for n ¼ 0: q 7.5. Distortion and dissipation Results (66) and (67) can be used to implement definitions (54) and (55) for the numerical normalized damped

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

1.0

ρ(A)

ξ = 0.00

0.8

N=2 4

1.0

ρ(A)

ξ = 0.04 NMK & OVIP

0.8

8

0.6

0.6

0.4

0.4

3667

N= 2 4

NMK, OVIP, Modal & Exact: ρ(A) = 1

0.2 0.0 -3.0

1.0

-1.8

-0.6

0.2

log(Ω) 0.6

1.8

3.0

ξ = 0.10

ρ(A)

8

0.0 -3.0

1.0

0.8

NMK & OVIP

0.8

0.6

Exact & N = 4, 6, 8

0.6

N=2

0.4

log(Ω) -1.8

-0.6

0.6

-0.6

0.6

1.8

ρ(A)

log(Ω) 3.0

ξ = 0.40 OVIP NMK

0.4

0.2 0.0 -3.0

-1.8

1.8

3.0

0.2

Exact & Modal

0.0 -3.0

-1.8

2 4

log(Ω) -0.6

0.6

1.8

3.0

Fig. 1. Estimates for the numerical spectral radius.

frequency and for the damping ratio and obtain the following expressions:  d ¼ Xd þ OðX2pþ2 Þ; X  n ¼ n þ OðX2pþ2 Þ:

ð69Þ ð70Þ

The values obtained for the distortion and dissipation measures defined by Eqs. (58) and (59) are shown in Figs. 2 and 3, respectively. It is noted that the values obtained with high-degree bases are affected by numerical round-off (in double precision). Moreover, the instability shown in Figs. 2 and 3 for higher frequencies is due to quadrant-change

0

ξ = 0.00

log(εΩ)

sensitivity in the calculation of the normalized numerical angular frequency (log(p)  0.5). This source of instability is weaker in the results presented in Fig. 3, which show that the numerical damping tends to zero as the frequency increases. One of the features enhanced by the tests reported above is the level of stability and accuracy that can be attained with the time integration method using relatively large time steps, that is, for relatively large values of X = xDt. Results (68)–(70) justify the quality of the values obtained for relatively small frequencies using low-degree bases. According to the study reported in Ref. [11], results (69) and (70) with

0

ξ = 0.04

log(εΩ)

NMK

-4

-4 NMK OVIP

2

-8

-8

2 8

8

-12

OVIP & Modal recover Exact

6

-12

N=4

-16 -3.0

0

-1.8

log(Ω) -0.6

0.6

1.8

3.0

ξ = 0.10

log(εΩ)

-16 -3.0

0

-4

6 N=4

-4 2

OVIP

-8

NMK OVIP

-12

6 N=4

-16 -3.0

-1.8

0.6

1.8

3.0

ξ = 0.40

2 8

8

-12

-0.6

log(εΩ)

NMK -8

log(Ω)

-1.8

log(Ω) -0.6

0.6

1.8

3.0

6 N=4

-16 -3.0

-1.8

log(Ω) -0.6

0.6

 d =Xd . Fig. 2. Solutions for the dispersion measure eX ¼ 1  X

1.8

3.0

3668

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

p = 0 define the numerical estimates for the damped frequency and for the damping ratio obtained with the Newmark (trapezoidal) rule, which is associated with a spurious (non-null, real, sign-unrestricted) eigenvalue, k1 6¼ 0.

T n ðtÞvn þ vp ðtÞ;

tion matrix, defined by Eq. (61), is recovered with the approximations corrected with the initial solution. The initial solution correction has no effect, also, on the estimates presented in Fig. 2 for the numerical spectral radius and on the results presented in Figs. 3 and 4, on the dispersion and dissipation measures. However, the modal particular solution recovers the exact values for the spectral radius, the damped frequency and the damping ratio, independently of the dimension of the polynomial time approximation basis. This should be expected, as the exact solution is now contained in finite element approximation. It is recalled that the exact values for the stability measures are recovered, also, using the generalized virtual impulse (GVIP) method, as it is based on the exact solution of the homogeneous scalar problem.

T n ðtÞan þ ap ðtÞ:

7.7. Convergence

7.6. Effect of corrective terms Approximations (32)–(34) are extended as follows to assess the influence of corrective terms in the stability and convergence of the solution: uðtÞ ¼

N X

T n ðtÞun þ up ðtÞ;

n¼1

vðtÞ ¼

N X n¼1

aðtÞ ¼

N X n¼1

These corrective terms are a quadratic solution that fits the initial conditions of the problem, and the mode decomposition solution (see Section 4 and Eq. (46), where now an = Dt1 for simplicity):

Convergence under both p- and h-refinement procedures is illustrated using a scalar second-order problem (1)–(3) with a Legendre polynomial solution of degree k = 7:

1 up ðtÞ ¼ u0 þ tv0 þ t2 a0 ; 2 N sys X up ðtÞ ¼ exp ððnn xn þ jxdn ÞtÞ/n uPn :

In the present context, p-refinement is implemented by increasing the degree p of the polynomial time approximation basis (with dimension N = p + 1). In the solutions shown in Fig. 4, the dimension of the basis is increased from N = 2, a linear approximation, to N = 8, the seventh degree approximation that recovers the exact solution (71). The alternative h-refinement procedure consists in decreasing the time step, Dt. The results shown in Fig. 4a are obtained using a single time step, Dt = 2, and the results presented in Fig. 4b are determined implementing an incremental process on four time steps, Dt = 0.5.

n¼1

When the approximations thus extended are implemented on the scalar problem (47), it is found that they do not affect definition (62) for the components of the forcing load vector. The exact amplification matrix, defined in Appendix 1, is recovered for the approximations corrected by the modal decomposition, and the approximate amplifica-

0

ξ = 0.01

log(εξ)

0

NMK OVIP

-4

uðtÞ ¼ P k ðtÞ;

-4

8 -8 -12

-4

-8

N=4

-12

log(Ω)

-16 -3.0

0

6

-1.8

-0.6

0.6

3.0

ξ = 0.10

log(εξ) NMK OVI P

1.8

2 8 6

N=4

-1.8

N=4

log(Ω) 0.6

1.8

3.0

ξ = 0.40

NM K OVI P

2 8

6

6

-12

log(Ω) -1.8

-0.6

log(εξ)

-8

-12 -16 -3.0

NM K OVI P

-4

-8

ð71Þ

ξ = 0.04

log(εξ)

-16 -3.0

0

2 8

1 6 t 6 þ1:

-0.6

0.6

1.8

3.0

-16 -3.0

N=4

log(Ω) -1.8

-0.6

0.6

Fig. 3. Solutions for the dissipation measure en ¼ 1  n=n.

1.8

3.0

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

a(t)/a(1)

a(1) = 378 1.0

0.5

0.5

0.0

0.0 Exact & N = 8 N=4

-0.5

N=6 N=2

-0.5

t

-1.0 -1.0

-0.5

0.0

0.5

v(t)/v(1)

1.0

v(1) = 28

0.5

0.5

0.0

0.0

-0.5

-0.5

t 0.0

0.5

1.0

t -0.5

0.0

0.5

u(t)

1.0

v(1) = 28

Exact

N=2

N=4

t

-1.0 -1.0

-0.5

0.0

0.5

1.0

u(t)

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0 -1.0

N=2

v(t)/v(1) 1.0

-0.5

Exact

-1.0 -1.0

1.0

-1.0 -1.0

a(1) = 378

a(t)/a(1)

1.0

3669

t -0.5

0.0

0.5

1.0

(a) Single-step solution, Δt = 2.0

-1.0 -1.0

Exact N=4

N=2 N=6

t -0.5

0.0

0.5

1.0

(b) Four-step solution, Δt = 0.5

Fig. 4. Solution estimates for the Legendre test (71) using polynomial bases of degree N  1.

The quality of the solutions thus obtained is insensitive to the relative values used for the mass, damping and stiffness coefficients, which define the amplitude of the forcing load that equilibrates the elastic, damping and inertia forces due to the displacement field (71) and to the associated (sixth-degree) velocity and (fifth-degree) acceleration fields. It is useful to analyse the estimates obtained at the endpoints of each time step and within each time interval. The graphs confirm that the approximate solutions do not enforce the initial condition set in each time increment, see Eq. (24). However, the error converges to zero as the (h- or p-) quality of the approximation is improved. They show, also, that convergence to the exact solution at the end of each time increment is monotonic and relatively faster. Moreover, and as it should be expected, the exact solution at the end-points is recovered as soon as the (polynomial) time basis accommodates the degree of the solution, that is, when N = 6 for the acceleration estimate and N = 8 for the velocity and displacement estimates. The same reason justifies the fact that the exact solutions within the time interval are recovered only when the time

basis accommodates the appropriate degree of approximation. This is clearly illustrated for the single-step solution, which shows, also, that the estimates obtained within the time interval can be rather poor. This is expected, as a polynomial basis of a given degree cannot recover a solution of higher degree defined on the same interval. It is stressed, also, that the relative rates of convergence found for the acceleration (highest rate), velocity and displacement (lowest rate) cannot be extrapolated. They are a direct consequence of the particular nature of this test, in which a polynomial basis, with the same degree in the approximation for the displacement, velocity and acceleration, see Eqs. (4), (8) and (9), is used to recover the exact solution of a polynomial problem. In summary, the results presented in Fig. 4 show that the solution is stable when implemented on large time steps and confirm that the estimates obtained at the end of the step converge to the exact solution, an essential feature for step-by-step integration methods. They indicate also that a rich approximation (e.g. wavelet) basis must be used to obtain adequate estimates of the exact solution within the time increment, when a large time step strategy is used.

3670

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678



7.8. Numerical tests The first test taken from Ref. [14] (with the correction of a typing mistake in the definition of the forcing load vector), models the forced vibration of an undamped spring– mass system, initially at rest and with a strongly distorted stiffness: 2

1 6 40

0 1

0

0

38 9 2 6 0 > = < a1 > 6 7 0 5 a2 þ 4 2 > ; : > 2 0 a3

38 9 8 9 2 > = = > < 10 > < u1 > 7 2 5 u2 ¼ 0 : > ; ; > : > : > 2 204 0 u3 2 6

The second test problem taken from the same reference models a strongly non-proportionally damped system, with initial conditions uT0 ¼ f 0:1 0:2 g and vT0 ¼ f 0:3 0:5 g:

10.0

2

2

3 ¼

2.5



 þ

4:24

3:70



3:70 4:35 a2

10 sinð0:702t  p=2Þ : sinð1:92606t þ p=2Þ

v1



v2

 þ

6

2

2

6



u1



u2

a1

"Exact" solutions: Exact & present (T=0.25, N=2)

1.5

2.0

0.5

-2.0

-0.5

-6.0

-1.5

-10.0 0.0 5.0

a1

The results presented in Fig. 5 compare the exact solution of the distorted stiffness problem with the estimates obtained with the Newmark (NMK) method and with the procedure described above implemented on a linear time basis, N = 2, for the same small step Dt = T = 0.2. This test is also used to illustrate the convergence of the solutions obtained by increasing the degree of the approximation to cubic, N = 4, and quintic bases, N = 6, under a relatively large time integration step, Dt = 2.0. The estimates obtained for the solution of the non-proportionally damped problem are shown in Fig. 6. They are obtained with linear and cubic bases for two time steps,

a1

6.0



5

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

18.0

t

-2.5

20.0

-3.5

t 0.0

v1 1.5

3.0

2.5

5.0

v1

7.5

10.0

12.5

15.0

17.5

20.0

22.5

25.0

"Exact" solutions: Exact & present (T=0.25, N=2)

1.0 1.0 0.5 -1.0 0.0 -3.0

-0.5

t

-5.0 0.0

4.0

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

18.0

-1.0

20.0

t

-1.5

u1

0.0

2.5

5.0

7.5

10.0

12.5

15.0

17.5

20.0

22.5

25.0

u1 "Exact" solutions: Exact, present (T=0.25, N=2) and OVIP (T=0.25) 1.5

3.2

1.0

2.4

0.5 1.6 0.0 0.8

-0.5

t

0.0 0.0

2.0

4.0

6.0

8.0

10.0

12.0

14.0

16.0

18.0

-1.0

20.0

t

-1.5 Exact & present (T=0.2, N=2)

NMK (T=0.2)

Present (T=2.0, N=6) Present (T=2.0, N=2)

Present (T=2.0, N=4) NMK (T=2.0)

Fig. 5. Estimates of the response of the spring–mass system with distorted stiffness.

0.0

2.5

5.0

7.5

"Exact" solutions Present (T=1.5, N=2)

10.0

12.5

15.0

17.5

NMK (T=0.25) Present (T=1.5, N=4)

20.0

22.5

25.0

NMK (T=1.5) GVIP (T=1.5)

Fig. 6. Estimates of the response of the non-proportionally damped forced system.

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

Dt = 0.25 and 1.5, and compared with the estimates obtained with the Newmark method and with the original and generalized versions of the virtual impulse method, denoted by OVIP and GVIP, respectively. The softening and hardening spring test problems used in Ref. [14], €u þ 100 tanhðuÞ ¼ 0 2

€ u þ ð1 þ u Þu ¼ 0

_ with uð0Þ ¼ 0 and uð0Þ ¼ 10;

_ with uð0Þ ¼ 0 and uð0Þ ¼ 60

are adopted here to illustrate the application of the time integration procedure to the modelling of nonlinear responses under intentionally unfavourable implementation conditions, particularly in what concerns the second test. A pedestrian method is used in detriment of gradient-improved convergence techniques, as the issue is the robustness of the time integration scheme in terms of convergence. The iteration hinges fully on the resulting non-residual term as the initial stiffness of the system is used throughout the process €unþ1 þ K 0 unþ1 ¼ F NLn to yield K0 = 100 and F NLn ¼ 100un  100 tanhðun Þ for the softening spring test, and K0 = 1 and F NLn ¼ u3n for the hardening spring test. Moreover, a time approximation basis best suited for large time stepping is used. The wavelet basis defined in Section 8 is applied, instead of a low-degree polynomial basis, which is more suitable to accommodate the (small) time increments adequate to the implementation of Newton-type techniques, Ref. [3]. The results shown in Fig. 7 are obtained setting a tolerance of 106 on the relative error on each order of the displacement approximation (4) and using time steps Dt = 0.75 and Dt = 0.15 for the softening and hardening spring models, respectively. Solutions found with the Newmark, HHT and VIP integration schemes, with increment Dt = 0.01, are reported in Ref. [14]. 8. Hyperbolic and parabolic problems The two alternative approaches in the solution of hyperbolic/parabolic problems based on the separation of variu

"Exact"

1.5

ables in space and time are briefly recalled next. The first approach, and the most commonly used, of relying on a problem-independent approximation in space, is recalled with the simple purpose of clarifying the application of the time integration procedure discussed here. The second approach, typical of spectral analysis, consists in discretizing the problem first in the time domain to obtain an elliptic Helmholtz-type problem, which is then solved using the procedures for approximation in space that are available. This is the approach that has motivated the development of the time integration procedure reported here, to support the solution of hyperbolic and parabolic problems using the Trefftz variant of the finite element method. Its implementation is illustrated using two testing problems. 8.1. Notation The equations governing an hyperbolic problem, defined in domain V with boundary S = Sr [ Su, and referred to a Cartesian system x and a time frame t, can be stated as follows: _ tÞ L kLdðx; tÞ þ fðx; tÞ ¼ m€dðx; tÞ þ cdðx; T N kLdðx; tÞ ¼ tðx; tÞ on S r ;  tÞ on S u ; dðx; tÞ ¼ dðx;

0.3

2.0

-0.3

-2.0

-0.9

-6.0

t 0.5

1.0

1.5

2.0

2.5

3.0

(a) Softening spring model (Δt = 0.75).

"Exact"

Present

10.0 6.0

0.0

ð72Þ ð73Þ ð74Þ

In the terminology of solid mechanics, in the equation of motion (72) d is the displacement field, L is the differential compatibility (gradient) operator, L* is the differential equilibrium (divergence) operator (the conjugate of L in geometrically linear problems), f is the body force vector, k, m and c are the local stiffness, mass and damping matrices, and d_ and €d are the velocity and acceleration vectors. In the boundary equilibrium (Neumann) condition (73), matrix N collects the components of the unit outward normal and t is the prescribed force vector, and d is the prescribed displacement vector in the boundary compatibility (Dirichlet) condition (74). The results presented below can be adapted to

0.9

-1.5

in V ;

dðx; 0Þ ¼ d0 ðxÞ in V ; _ 0Þ ¼ d_ 0 ðxÞ in V : dðx;

u

Present

3671

t

-10.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

(b) Hardening spring model (Δt = 0.15).

Fig. 7. Nonlinear free vibration of single-degree-of-freedom systems.

3672

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

Dirichlet conditions written in terms of prescribed velocities. 8.2. Problem-independent approximation bases For a typical (conform) displacement element, with domain Ve and boundary Se, the displacement, velocity and acceleration fields are approximated in form dðx; tÞ ¼ UðxÞuðtÞ in V e ; _ tÞ ¼ UðxÞvðtÞ in V e ; dðx; € tÞ ¼ UðxÞaðtÞ in V e dðx;

un ¼ un

_ with vðtÞ ¼ uðtÞ and a(t) = u¨(t) and where the space b approximation basis is real, in general, UðxÞ ¼ UðxÞ, and chosen to ensure that interelement and boundary displacement continuity conditions are locally satisfied, while letting: uðx; 0Þ ¼ UðxÞu0 in V e ; vðx; 0Þ ¼ UðxÞv0 in V e : Different techniques can be used to discretize the hyperbolic problem into the second-order problem (29), to yield the following expressions for the mass, damping and stiffness matrices and for the force vector: Z M ¼ UT mU dV e ; Z C ¼ UT cU dV e ; Z K ¼ ðLUÞT kðLUÞ dV e ; Z Z T e F ¼ U f dV þ UTt dS e : The solving system expression (40) holds with the load decomposition defined by Eq. (42), and the estimates for the time components of the displacement, velocity and acceleration fields are obtained using Eqs. (32), (43) and (44), respectively, according to the procedure summarized in Section 5.5. 8.3. Problem-dependent approximation bases In the extensions (75)–(77) of the basic time approximations (32)–(34), the generalized displacement, velocity and acceleration vectors, un(x), vn(x) and an(x), define the space-dependent weights of time approximation function Tn, with x 2 V: dðx; tÞ ¼

N X

Tn ðtÞun ðxÞ;

ð75Þ

Tn ðtÞvn ðxÞ;

ð76Þ

Tn ðtÞan ðxÞ:

ð77Þ

n¼1

_ tÞ ¼ dðx;

N X n¼1

€ dðx; tÞ ¼

N X n¼1

The velocity and acceleration integration rules (38) and (39) still hold for multi-dimensional problems in the space domain and the hyperbolic problem defined by Eqs. (72)– (74) is replaced by the following set of N uncoupled (Helmholtz-type) elliptic problems   Xn Xn  m þ c un  f 0n in V ; ð78Þ L kLun þ f n ¼ Dt Dt NT kLun ¼ tn on S r ð79Þ ð80Þ

on S u ;

where the body force vector and the prescribed force and displacement vectors are defined by Z Dt N X H 1 nm f n ðxÞ ¼ Tb m fðx; tÞdt; ð81Þ Dt 0 m¼1 Z Dt N X H 1 nm tn ðxÞ ¼ Tb mtðx; tÞdt; ð82Þ Dt 0 m¼1 Z Dt N X H 1 nm un ðxÞ ¼ Tb m uðx; tÞdt ð83Þ Dt 0 m¼1 and where the force vector associated with the initial conditions of the problem is defined by   X0n Xn X0n _ 0 m þ c d0 þ md0 : fn ¼ Dt Dt Dt The so-called Trefftz constraint consists in constraining the selection of the space approximation basis, say Un(x), un ðxÞ ¼ Un ðxÞqn

in V e

ð84Þ

to the set of formal solutions of the homogeneous governing differential equation (78):   Xn Xn  m þ c Un in V e : L kLUn ¼ Dt Dt The solution of problem (78)–(80) using this space discretization criterion can be found, for instance, in Ref. [2]. Two alternative models are developed, namely the displacement and the stress models, which are designed to produce (weakly) conform and equilibrated solutions, respectively. They combine features that typify the finite element method, namely domain decomposition, sparsity and symmetry, with the main feature of the boundary element method, as all coefficients present in the algebraic solving system are defined by boundary integral expressions. The added advantage of the Trefftz finite element approach is that regular approximation bases are used instead of singular (fundamental solution) bases that typify the implementation of the boundary element method. Two final aspects are worth mentioning. The first is that problem (78)–(80) reduces to the Helmholtz problem of frequency domain analysis for periodic problems (with null forces associated with the initial conditions, as noted in Section 4.3). The costs of using spectral techniques in the analysis of non-periodic problems by periodic extension are well known and motivated the non-periodic spectral

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

extension suggested here to support the implementation of the Trefftz method in the solution of time domain problems. The second aspect worth mentioning is that the implicit methods of time integration that have been proposed, as reviewed for instance in Ref. [14], can be used to discretize first the hyperbolic problem (72)–(74) to obtain the elliptic form (78)–(80) needed to implement the Trefftz variant of the finite element method, as first suggested by Jirousek and Qin [8]. This is equivalent to implement a spectral decomposition in time with a single frequency spectrum, N = 1 in approximations (75)–(77), which has a negative effect in the rate of convergence of the Trefftz finite element solution and forces the use of small time-stepping strategies. Still in the context of linear problems, the relative advantages offered by the integration method presented here are clearly demonstrated by the two linear applications reported next. 8.4. Bar subject to impact load The bar used in the undamped tests reported below has length L = 100 mm, modulus of elasticity E = 2  104 MPa and specific mass q = 2  104 N s m4. It is fixed at the origin and the end-section is subject to an impact load with intensity r0 = 2 MPa. The axial stress and displacement estimates, r and u, are normalized to the scaling values r0 and u0 = r0L/E, respectively, and the non-dimensional pffiffiffiffiffiffiffiffi ffi time parameter used is defined by: s ¼ E=q  t=L. The objective of the first set of tests is to compare the solutions obtained with the Newmark method and with the weakest approximation (linear, N = 2) of the time integration procedure described here, using in both cases linear displacement (constant stress) elements. The results presented in Fig. 8 show that the stress wave estimate obtained with the Newmark procedure remains relatively poor when a very small time increment is used,

3.0

τ = 0.25

σ/σ0

Ds = 0.0625, together with a high refinement of the bar in 64 (uniform) elements. These results show that the scheme reported here is only marginally affected by spurious local oscillations, which visibly pollute the Newmark estimates. Obviously, these estimates can be improved using alternative, and by now well-tuned, lower-order time-stepping algorithms, and, also, by improving the quality of the approximation in space. This form of refinement has a visible effect on the quality of the estimates produced by the procedure reported here. The last test on the bar subject to an impact load is designed to illustrate the ‘‘spectral” capabilities of the time integration procedure presented here, under no periodicity constraints. The same problem is solved in a single time step, Ds = 2.0, and a single Trefftz finite element, which is derived for the ‘‘forcing frequencies”, Xn, generated by a wavelet time approximation basis. The time basis used in the implementation of the Trefftz element is built on N = 32 scaling functions adapted to the (unit) interval, for family 4 and refinement 5. The results obtained are compared in Figs. 9 and 10 with the Fourier solution implemented on a periodicity window Ds = 4.0 and using N = 100 approximation modes. It is noted that the incremental linear displacement solutions are obtained by solving, at each time step, a system of Nsys = 64 equations, while the ‘‘spectral” solutions presented in Figs. 9 and 10 are obtained solving N/2 = 16 (complex) scalar systems, as the discretization into a single Trefftz element leads to a single degree-of-freedom system, the displacement at the free-end of the bar. This crude approximation in the space dimension and the relative weakness of the wavelet bases are responsible for the difficulty of modelling the propagation of the stress wave front, clearly shown in Fig. 10. Improved convergence is obtained by refining the approximation in space and time, as shown below.

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

τ = 0.50

σ/σ0

1.0 Exact

0.5 0.0 0.000

3.0

3673

Newmark

Present

0.5

x/L 0.125

0.250

0.375

0.500

0.625

0.750

0.875

1.000

τ = 0.75

σ/σ0

0.0 0.000

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5 0.0 0.000

x/L 0.125

0.250

0.375

0.500

0.625

0.750

0.875

1.000

τ = 1.00

σ/σ0

0.5

x/L 0.125

0.250

0.375

0.500

0.625

0.750

0.875

1.000

0.0 0.000

x/L 0.125

0.250

0.375

0.500

0.625

Fig. 8. Stress estimates for incremental analysis with Ds = 0.0625 and 64 elements.

0.750

0.875

1.000

3674

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

2.00 1.75

Exact and present solutions

u/u0

2.00

τ = 2.00 1.75

Linear approximation in time Linear displacement elements

Exact and present solutions

u/u0

1.75

τ = 2.00 1.75

Wavelet approximation in time Single Trefftz displacement element

1.50

1.50

1.50

1.25

1.25

1.25

1.25

1.00

1.00

1.00

1.00

0.75

0.75

0.75

0.75

0.50

0.50

0.50

0.50

0.25

0.25

0.25

0.25

1.50

x/L 0.00 0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

x/L 0.00 0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000

(a) Analysis with Δτ = 0.0625 and 64 elements.

(b) Single-time step and single-element analysis.

Fig. 9. Displacement estimates for the impact load test.

3.0

τ = 0.25

σ/σ0

3.0

2.5

2.5

2.0

2.0

1.5

1.5

Exact solution Trefftz element (32 modes)

1.0

1.0

Fourier solution (100 modes)

0.5 0.0 0.000

3.0

0.125

0.250

0.375

0.500

0.625

0.5 0.750

σ/σ0

0.875

1.000

t = 0.75

0.0 0.000

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.0 0.000

x/L 0.125

0.250

0.375

0.500

0.625

0.750

σ/σ0

0.875

1.000

0.250

0.375

0.500

0.625

0.750

σ/σ0

t = 1.25

0.0 0.000

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.875

1.000

t = 1.00

x/L 0.125

0.250

0.375

0.500

0.625

0.750

σ/σ0

0.875

1.000

t = 1.50

0.5

0.0 0.000

3.0

x/L 0.125

0.5

0.5

3.0

τ= 0.50

σ/σ0

x/L 0.125

0.250

0.375

0.500

0.625

0.750

σ/σ0

0.875

1.000

t = 1.75

0.0 0.000

3.0

2.5

2.5

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.5

0.0 0.000

x/L 0.125

0.250

0.375

0.500

0.625

0.750

0.875

1.000

x/L 0.125

0.250

0.375

0.500

0.625

0.750

σ/σ0

0.0 0.000

0.875

1.000

t= 2.00

x/L 0.125

0.250

0.375

0.500

0.625

0.750

0.875

1.000

Fig. 10. Stress estimates for the impact load test (single-time step, single-element mesh).

8.5. Stress relaxation of soft tissue sample The unconfined compression test on cartilage specimens represented in Fig. 11 is often used to illustrate the perfor-

mance of hybrid and mixed finite element formulations developed to model the response of incompressible biphasic media. The data used to implement this two-dimensional parabolic problem, m = 0 in Eq. (72), is taken

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678 y

h

Modelled quadrant

u x

to

tmax

w

Fig. 11. Unconfined compression (plane strain) test on incompressible soft tissue.

from Vermilyea and Spilker [17]: width w = 6.35 mm, height h = 1.78 mm, modulus of elasticity E = 0.675 MPa, Poisson’s ratio t = 0.125, permeability j = 7.6  1015 m4 N1 s1, and fluid fraction /f = 0.83. The displacement-driven, ramp-load programme shown in the same figure is designed to reach a deformation of 5% at instant t0 = 500 s, which corresponds to a prescribed displacement  u ¼ 0:0445 mm. The duration of the loading is tmax = 1000 s and two boundary conditions are tested, for adhesive and lubricated loading platens, respectively. The time integration method reported here is implemented on a wavelet basis with dimension N = 256 (family 3 and refinement 7). It is applied in a single time step, as its support is identified with the full duration of the test, Dt = tmax in Fig. 11. Hybrid-Trefftz stress elements are used to model a quadrant of the specimen. Two coarse meshes are used, namely

ν/νo

a single-element and a regular 2  2-element mesh. In the adhesive platen test, the corresponding solving systems have a total of 100 and 412 degrees-of-freedom, respectively. The numbers are marginally higher for the lubricated platen test, 106 and 424, respectively. It is noted that each of these dimensions combines domain (stress modes) and boundary (displacement modes) degrees-offreedom, as local condensation at element level is not called upon. Moreover, no use is made of the adaptive refinement capabilities of the approximation bases used, either in time (wavelet basis) or in space (Trefftz basis). The evolution in time of the results obtained at control point x = w/2 and y = 0, see Fig. 11, with the single-element mesh and with the regular 2  2-element mesh, in a single time step, Dt = 103 s, is shown in Fig. 12. The velocity components shown are in the dominant direction of the flow. The stress components are measured in the solid skeleton and p denotes the total pressure. The graphs show that the single-element solutions model captures adequately the peaks in acceleration occurring at the end-points of the loading phase. The adhesive platen test solution obtained with the single-element mesh recovers the values reported by Vermilyea and Spilker [17] using a trapezoidal rule with time step Dt = 5 s and a regular mesh of 12  6 pairs of hybrid elements. The animation of the displacements, in the fluid and in the solid phases, and of the stress and pressure components during the full loading process can be accessed using the address www.civil.ist.utl.pt/HybridTrefftz, under link ‘Animation of the response of hydrated soft tissues’.

vo = 0. 6 μ m/s & to= 500 s

σ/σ o

0.4

1.0 0.8

Velocity in fluid

σ o = 50 kPa & t o = 500 s σ xxs / σo , σ xys / σo , p /σ o ≈ 0

0.2

0.6

Control point: x = 0. 5w ; y = 0

0.0

0.4

Single-element mesh

-0.2

Source point: x = 0. 5w ; y = 0

-0.4

2x2 element mesh

Velocity in solid

0.1 -0.1

-0.6

-0.3

-0.8

-0.5 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75

3675

s σ yy / σo

-1.0

t / to

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75

t / to

(a) Unconfined (lubricated) compression test. 1.3

ν/νo

v o = 0. 6 μ m/s & to= 500 s

Velocity in fluid

0.4

Contr ol point: x = 0.5w ; y = 0

0.0

Single-element mesh

-0.2

0.1 -0.2

σ/σ o

0.2

1.0 0.7

0.4

Velocity in solid

-0.8 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75

σ xxs / σo , σ xys / σo , p /σ o ≈ 0

-0.4

Source point: x = 0.5w ; y = 0

-0.6

2x2 element mesh

-0.8

-0.5

σ o = 50 kPa & t o = 500 s

σ yys / σo

-1.0

0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 t / to t / to (b) Unconfined (adhesive) compression test.

Fig. 12. Velocity and stress estimates for the soft tissue modelling problem.

3676

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

9. Closure The Galerkin, finite element time integration procedure reported here can be implemented in both incremental and spectral forms, in the latter case released from the periodicity (or extended periodicity) conditions that constrain Fourier analysis. It is open to alternative time bases, independently of the approximation used in the discretization in space. It is unconditionally stable, insensitive to overshooting and may yield competitive convergence rates, either by reducing the time step or by enriching the approximation basis. This procedure replaces the solution of a system with dimension Nsys by the solution of an equivalent set of fully uncoupled N algebraic systems with the same dimension, Nsys, with N representing the dimension of the time basis. This dimension is independent and usually much smaller than the dimension of the system being solved, and does not call upon its modal decomposition. The uncoupling condition is attained solving a supporting eigenvalue problem, with the dimension of the time approximation basis, which generates pairs of complex conjugate Fourier-equivalent spectral frequencies, while preserving the properties of the (not necessarily symmetric/Hermitian) system matrix. This procedure cannot compete with the most commonly used lower-order implicit methods for small time stepping, but offers advantages when the use of large time steps may be advisable, even under nonlinear conditions, see Ref. [3]. The implementation of the spectral decomposition of the time basis, which in the present context is numerically equivalent to solve N time-step problems, depends strongly on the richness (dimension) of the time approximation bases and on its sensitivity to ill-conditioning. Preliminary tests have shown that systems of wavelets accommodate these basic requirements for non-periodic spectral-type solutions. Parabolic/hyperbolic problems can be solved using the time integration procedure presented here, either by discretizing in the space domain first, as it is the common practice to obtain first-/second-order time-dependent systems, or by discretizing in space after discretizing in time the state equations. In this latter case, the parabolic/hyperbolic problem is replaced, in the present case, by N/2 pairs of complex conjugate Helmholtz-type elliptic problems, with wavenumbers that are dictated by the equivalent forcing frequencies generated by the time approximation basis. The end result is the same when the space approximation basis is problem-independent. However, when this dependency is exploited in the second approach, it becomes possible to extend into the time domain the advantages offered by the Trefftz variant of the finite element method in the solution of elliptic problems.

presentation of the paper. The author wishes, also, to thank Mr. I. Moldovan and Mr. M. Toma for preparing the nonlinear testing problems and the numerical test on soft tissue modelling, respectively. This research has been supported by Fundacßa˜o para a Cieˆncia e a Tecnologia. Appendix 1. Exact solution of scalar second-order problems Dropping the subscript, for simplicity of the notation, the solution of Eq. (45) is uðtÞ ¼ expðnxtÞðcosðxd tÞX 1 þ sinðxd tÞX 2 Þ þ up ðtÞ; Z t 1 up ðtÞ ¼ ðMxd Þ expðnxðt  sÞÞ sinðxd ðt  sÞÞf ðsÞds; 0

pffiffiffiffiffiffiffiffiffiffiffiffiffi where xd ¼ 1  n2 x is the damped frequency. The initial conditions yield the identifications: X 1 ¼ u0 X 2 ¼ x1 d ðnXu0 þ v0 Þ: Appendix 2. Exact amplification matrix and forcing load vector The solution defined in Appendix 1, with the notation introduced in Section 7 for the normalized angular frequency, X, and for the normalized damped frequency, Xd, can be used to derive the following expression for the amplification matrix present in Eq. (48), where definitions (50) and (51) are used: 2 3 Xd a þ nXb b 0 6 7 A ¼ X1 X2 b aXd  nXb 05 d 4 X2 Xd a þ nX3 b 2nXXd a  ð1  2n2 ÞX2 b 0

ð1:2Þ The expression for the forcing load vector, present in the same equation is X2 f L2 L1 L0 g; K Z 1 L2 ¼ sðsÞeðsÞf ðsÞds; LT ¼

L1 ¼

Z

0 1

ðcðsÞ  nXsðsÞÞeðsÞf ðsÞds; Z 1 ð2nXcðsÞ  ðX2d  n2 X2 ÞsðsÞÞeðsÞf ðsÞds; L0 ¼ f ðDtÞ  0

0

where s = t/Dt; e(s) = exp(nX(1  s)); (1  s)); s(s) = sin(Xd(1  s))/Xd.

c(s) = cos(Xd

Acknowledgements

Appendix 3. Properties of complete polynomial bases

The author wishes to thank the reviewers, as their comments were essential to improve both the structure and the

It is convenient, for simplicity, to construct the polynomial basis (see Section 4.2) on a complete set of orthonor-

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

mal Legendre polynomials, Pn(f), to yield the following results: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi un ðsÞ ¼ 2n  1P n1 ðfÞ with f ¼ 2s  1; Z 1 um un ds ¼ dmn ; 0 Z 1 u0m un ds ¼  dmn ; 0

where dmn is the Kronecker symbol, dmn ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð2m  1Þð2n  1Þ if m > n and m + n is odd, and  dmn ¼ 0 otherwise. In particular Z 1 um ds ¼ dm1 ; 0

Z

1

u0m ds ¼ um ð1Þ  um ð0Þ ¼  dm1 :

0

The procedure summarized in Section 4 can be used to verify that the following relations hold for a complete polynomial basis: N X

a0n b0n ¼ 1;

is obtained using the result d0m ¼ 1=m!, valid for a complete polynomial basis of degree p and m 6 2p + 1. Appendix 4. Power series approximation of amplification matrices The power series expansion of the exact amplification matrix, defined in Appendix 2 under notation (50) and (51), is expressed as follows when result (3.1) is used (seeAppendix 3): 1 X d n Xn An ; ðn þ 1Þ! n¼0 2 n þ 1 þ 2nX 1 6 2 An ¼ 4 X nþ1 2 2 ðn þ 1ÞX X  2nXðn þ 1Þ

A¼

Aa ¼

^ a0n X0n ¼ 1;

n¼1 ^ a0n Xn

2pþ1 X n¼0

X0n ¼ Xn b0n ; bkn Xn ¼ bkþ1 n ; The power series expansion of discriminant Dm, defined by Eq. (60) is 1 X n Dm ¼ d n ðX=Xm Þ d 0 ¼ 1; d 1 ¼ 2n; n¼0

d 2 ¼ ð1  4n2 Þ;

d 3 ¼ 4nð1  2n2 Þ;

d 4 ¼ 1  12n2 þ 16n4 ; 2

d 6 ¼ ð1  24n2 þ 80n4  64n6 Þ; . . .

ð3:1Þ

to yield the following expression for parameter ckm , defined by Eq. (63): 1 X ckm ¼ d n dkmþn Xn ; n¼0

T n ðDtÞbkn =Xmn :

The power series expansion of the exact solution for the components of the forcing load vector, defined in Appendix 2, is written in form 1 F X2 X d n Xn ; K n¼0 ðm þ n þ kÞ!

when result (3.1) is used and a polynomial load of degree k is assumed: k

F ðtÞ ¼ F ðt=DtÞ =k!: When this load is inserted in definition (62) and results (63)–(65) and (3.1) are used, the following expression is found for the approximate, polynomial solution:

n¼1

The expression that supports the accuracy and stability analyses, ckm ¼ c0mþk ¼

n¼0

d n Xn þ OðXM Þ ðm þ k þ nÞ!

with M ¼ 2ðp þ 1Þ  ðm þ kÞ

ð3:2Þ

OðX2pþ3 Þ 0

Appendix 5. Power series approximation of forcing load vectors

Lm ¼

4

d 5 ¼ 2nð3  16n þ 16n Þ;

M 1 X

7 0 5: 0

d n Xn  An þ A; ðn þ 1Þ!

OðX2pþ4 Þ

ckn ¼ ck1 nþ1 :

N X

3

where the residual term is defined by 2 3 OðX2pþ2 Þ OðX2pþ1 Þ 0 7  ¼6 A 4 OðX2pþ3 Þ OðX2pþ2 Þ 0 5:

¼ T n ðDtÞ;

dkm ¼ dmk ¼

0

The power series expansion of the approximate solution amplification matrix, defined by Eq. (61), is obtained using result (3.1), to yield,

n¼1 N X

3677

Lam ¼

1 X F X2 M d n Xn F X2 þ OðXM Þ K n¼0 ðm þ n þ kÞ! K

with M ¼ 2ðp þ 1Þ  ðm þ kÞ: The highest-degree term in a polynomial description of the forcing load is the only relevant term in the analysis of the truncation error of the approximate solution.

3678

J.A.T. de Freitas / Comput. Methods Appl. Mech. Engrg. 197 (2008) 3657–3678

References [1] I. Daubechies, Orthonormal bases of compactly supported wavelets, Commun. Pure Appl. Math. 41 (1988) 909–996. [2] J.A.T. Freitas, Hybrid-Trefftz displacement and stress elements for elastodynamic analysis in the frequency domain, CAMES 4 (1997) 345–368. [3] J.A.T. Freitas, Mixed finite element formulation for the solution of parabolic problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3425–3457. [4] J.A.T. Freitas, Time integration and the Trefftz Method: Part I – firstorder and parabolic problems; Part II – second-order and hyperbolic problems, CAMES 10 (2003) 453–477. [5] J.A.T. Freitas, J.P.M. Almeida, E.M.B.R. Pereira, Non-conventional formulations for the finite element method, Comput. Mech. 23 (1999) 488–501. [6] C. Hoff, P.J. Pahl, Development of an implicit method with numerical dissipation from a generalized single-step algorithm for structural dynamics, Int. J. Numer. Methods Engrg. 67 (1988) 367– 385. [7] J. Jirousek, Basis for development of large finite elements locally satisfying all field equations, Comput. Methods Appl. Mech. Engrg. 14 (1978) 65–92. [8] J. Jirousek, Q.H. Qin, Application of hybrid-Trefftz element approach to transient heat conduction analysis, Comput. Struct. 58 (1996) 195– 201.

[9] P. Monasse, V. Perrier, Orthonormal wavelet bases adapted for partial differential equations with boundary conditions, SIAM J Numer. Anal. 29 (1998) 1040–1065. [10] N.M. Newmark, A method of computation for structural dynamics, in: Proc. ASCE 85 EM3, 1959, pp. 67–94. [11] P. Pegon, Alternative characterization of time integration schemes, Comput. Methods Appl. Mech. Engrg. 190 (2001) 2707–2727. [12] J.P.P. Pina, J.A.T. Freitas, L.M.S. Castro, Use of wavelets in dynamic analysis, in: Me´todos Computacionais em Engenharia, APMTAC, Lisbon, 2004 (in Portuguese). [13] W. Ritz, Gesammelte Werke, Gauthier-Villars, Paris, 1911. [14] K.K. Tamma, X. Zhou, D. Sha, The time dimension: A theory towards the evolution, classification, characterization and design of computational algorithms for transient/dynamic applications, Arch Comput. Methods Engrg. 7 (2000) 67–290. [15] E. Trefftz, Ein Gegenstu¨ck zum Ritzschen Verfahren, in: Proc. 2nd Int. Cong. Appl. Mech., Zurich, 1926. [16] Trefftz Net. . [17] M.E. Vermilyea, R.L. Spilker, A hybrid finite element formulation of the linear biphasic equations for soft hydrated tissues, Int. J. Numer. Methods Engrg. 33 (1992) 567–594. [18] W.L. Wood, A unified set of single-step algorithms, Part 2: theory, Int. J. Numer. Methods Engrg. 20 (1984) 2303–2309. [19] O.C. Zienkiewicz, W.L. Wood, N.W. Hine, R.L. Taylor, A unified set of single-step algorithms, Part 1: General formulation and applications, Int. J. Numer. Methods Engrg. 20 (1984) 1529–1552.