A simplified variant of the time finite element methods based on the shape functions of an axial finite bar

Nội dung text: A simplified variant of the time finite element methods based on the shape functions of an axial finite bar

1. Journal of Science and Technology in Civil Engineering, HUCE (NUCE), 2021, 15 (4): 42–53 A SIMPLIFIED VARIANT OF THE TIME FINITE ELEMENT METHODS BASED ON THE SHAPE FUNCTIONS OF AN AXIAL FINITE BAR Thanh Xuan Nguyena,∗, Long Tuan Tranb aFaculty of Building and Industrial Construction, Hanoi University of Civil Engineering, 55 Giai Phong road, Hai Ba Trung district, Hanoi, Vietnam bVietnamese - German Construction Vocational Training Centre, College of Urban Works Construction, Model House A13, Yen Thuong street, Gia Lam district, Hanoi, Vietnam Article history: Received 02/8/2021, Revised 04/10/2021, Accepted 10/10/2021 Abstract In the field of structural dynamics, the structural responses in the time domain are of major concern. There already exist many methods proposed previously including widely used direct time integration methods such as ones in the β-Newmark family, Houbolt’s method, and Runge–Kutta method. The time finite element methods (TFEM) that followed the well-posed variational statement for structural dynamics are found to bring about a superior accuracy even with large time steps (element sizes), when compared with the results from methods mentioned above. Some high-order time finite elements were derived with the procedure analogous to the conventional finite element methods. In the formulation of these time finite elements, the shape functions are like the ones for a (spatial) 2-order finite beam. In this article, a simplified variant for the TFEM is proposed where the shape functions similar to the ones for a (spatial) axial bar are used. The accuracy in the obtained results of some numerical examples is found to be comparable with the accuracy in the previous results. Keywords: variational formulation; finite element; time finite element; non-linear axial bar element; shape function; structural dynamics; accuracy. © 2021 Hanoi University of Civil Engineering (HUCE) 1. Introduction The analysis of vibration response is of crucial importance in practical design of structures sub- jected to external time-varying actions. For obtaining of dynamic responses of a structure, the finite element method or a modal superposition approach is commonly used to spatially discretize the struc- ture, thus leading to a set of ordinary differential equations (ODE) in time whose primary unknowns are normally the time-varying nodal displacements of the structure. Practically, this set of ODEs can be solved by one of many time stepping approaches [1,2]. In general, there are mainly two classes among the direct time integration methods: explicit and implicit. Although requiring much more com- putations within a time step than that needed for an explicit method, implicit methods (such as ones in β-Newmark family, Houbolt’s method, and Runge-Kutta method) often possess unconditional sta- bility. To the contrary, explicit methods are conditionally stable. When coupled with the conventional ∗Corresponding author. E-mail address: thanhnx@nuce.edu.vn (Nguyen, T. X.) 42
2. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering finite element (FE) computation, the step size in any explicit method depends on the FE spatial mesh size and thus requires more computational effort in total for certain additional steps. For both implicit and explicit methods, a priori error analysis is often not easily available, since they are all derived in the spirit of finite difference. A review of implicit and explicit methods can be found in [3–5]. A different approach for solving dynamic problems is to use of the time (or temporal) finite ele- ment method (TFEM). The core idea of this approach is to discretize the time domain of the dynamic problem into “time” elements connected at chosen instants of time called as nodes. This approach has several potential advantages since the formulation of time finite element (TFE) can be applicable to both energy equation, and directly to the equations of motion. Usually, the TFEM approximation yields an accuracy superior to that of more conventional time stepping schemes at same computa- tional cost [6]. Furthermore, the formulation is convenient for computer implementation. The pioneer approaches, based on Hamilton’s Law of varying action, can date back to the research of Argyris and Scharpf, who employed Hermite cubic interpolation polynomials (akin to the beam finite ele- ment) to express the response over each time finite element [7]. Simkins [8] introduced all boundary and essential conditions into the “variational statement” as natural boundary conditions. This makes the variational statement suitable for obtaining approximate solutions for initial and boundary value problems. Employed this variational statement Simkins developed early finite element in time. Fried [9] applied the approach based on Hamilton’s principle to study the transient response of a damped system and transient heat conduction in a slab. Fried used a step-by-step approach to avoid storing and working with large matrices. Zienkiewicz and Parekh [10] used a time finite element approach to solve heat conduction problems. The formulation was based on Galerkin procedure over a time interval. In [11], Hulbert also employed the time-discontinuous Galerkin method and incorporates stabilizing terms having least-squares form. A general convergence theorem can be proved in a norm stronger than the energy norm. French and Peterson [12] proposed a time-continuous finite element method by transforming the second-order differential equations into first-order ones. Some other re- searchers have presented the variational formulation by allowing the TFE solution to be discontinuous at the end of each time element interval. Tang and Sun [13] introduced a unified TFE framework for the numerical discretization of ordinary differential equations based on TFE methods. In [6], Park extended a bi-linear formulation for linear systems and developed a TFEM to obtain transient re- sponses of both linear, nonlinear, damped and undamped systems. The sensitivity of the response with respect to various design parameters was also established. Through numerical examples, Park illustrated many applications of TFEM in sensitivity analysis, in parametric identification of non- linear structural dynamic systems, and in optimal control. Results for both the transient response and its sensitivity to system parameters, when compared to a previously available approach that em- ploys a multi-step method, are excellent. Recently, Wang and Zhong [14] proposed a time continuous Galerkin finite element method for structural dynamics. Its convergence property was proved through an a priori error analysis. The method provided by Wang and Zhong can give very accurate results compared with many other earlier methods. In 2020, Nguyen and Tran [15] proposed several devel- opments of high-order TFEs including the bp-TFEs which are analogous to the spatial second-order beam element, with the time-to-go (T − t) raised to the power of p. The formulation of bp-TFEs is based on the six shape functions of polynomials of degree 5 leading to element matrices of big sizes. Thus, in this article, we proposed a simplified variant of TFEM that can reduce the size of element matrices, using four shape functions for an axial bar element which are polynomial of degree 3 only. The accuracy in the results obtained by using this type of TFEs is less than that by using bp-TFEs, but is still superior than that by methods in β-Newmark family and is comparable to the ones given by 43
3. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering [14] in many cases. The rest of the article is as follows. In Section 2, a brief introduction of the variational formulation and bp-TFEs is presented. In Section 3, we propose a simplified variant of high-order TFEs, called as ap-TFE where “a” reminds the analogy to “axial” bar, whose shape functions of axial finite bars are used. Numerical examples are shown in Section 4, to illustrate the use of proposed ap-TFE. The article is finalized with conclusions in Section 5. 2. Variational formulation and bp-TFEs [15] Consider a structure having n degrees of freedom in the time interval LT =]0, T[. The mass n×n n×n matrix M ∈ R , and stiffness matrix K ∈ R are, generally, symmetric and positive definite. The n×n damping matrix C ∈ R is, generally, symmetric non-negative definite. The structure is subjected 2 2 to initial conditions u0 and u˙ 0, in addition to an external load F ∈ L (LT ), where L (LT ) is denoted for the Hilbert space on the time interval LT . The governing equation with the unknown u is as follows Lu ≡ Mu¨ + Cu˙ + Ku = F, u(0) = u0, u˙(0) = u˙ 0 (1) 2 Denote H (LT ) as the Sobolev space of order two. The two spaces were introduced as follows 2 n 2 o H0p (LT ) = u ∈ H (LT ) : u(0) = u0, u˙(0) = u˙ 0 (2) 2 n 2 o H00 (LT ) = u ∈ H (LT ) : u(0) = 0, u˙(0) = 0 (3) The strong form of structural dynamics in the above equation is equivalent to the following formula- 2 tion [14]: find u ∈ H0p(LT ) such that 2 B(u, v) = `(v), ∀v ∈ H00 (LT ) (4) where Z T Z s Z T B(u, v) = v˙T (Mu¨ + Cu˙ + Ku)dtds = (T − t)v˙T (Mu¨ + Cu˙ + Ku)dt (5) 0 0 0 Z T Z s Z T `(v) = v˙T Fdtds = (T − t)v˙T Fdt (6) 0 0 0 Let the interval [0, T] be divided into a finite number N of non-overlap sub-intervals each of which 2 2 2 has the length of Ti. Denote H0p,τ(LT ) and H00,τ(LT ) as the finite dimensional sub-spaces of H0p(LT ) 2 2 2 and H00(LT ), respectively. Usually, H0p,τ(LT ) and H00,τ(LT ) are assumed to be of polynomial forms in each element with degree greater than or equal to two. Then, by referring to the variational formulation 2 (5), the problem can be stated as: find uτ ∈ H0p,τ(LT ) such that 2 B (uτ, vτ) = ` (vτ) , ∀vτ ∈ H00,τ (LT ) (7) Based on the above variational formulation, the work in [15] proposed high-order bp-TFEs whose element stiffness matrix K and nodal element load vector f are given as follows Z Ti   K = (T − t)p H˙ T MH¨ + CH˙ + KH dt (8) 0 Z Ti f = (T − t)p H˙ T Fdt (9) 0 44
4. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering where a family of bp-TFEs can be obtained by selecting different values of p. The matrix H is analo- gous to the matrix of shape functions of a compatible second-order beam element with three equidis- tant nodes. 2 3 4 5  2 3 4 5 H1 = 1 − 23ξ + 66ξ − 68ξ + 24ξ , H2 = ξ − 6ξ + 13ξ − 12ξ + 4ξ Ti 2 3 4  2 3 4 5 H3 = 16ξ − 32ξ + 16ξ , H4 = −8ξ + 32ξ − 40ξ + 16ξ Ti (10) 2 3 4 5  2 3 4 5 H5 = 7ξ − 34ξ + 52ξ − 24ξ , H6 = −ξ + 5ξ − 8ξ + 4ξ Ti where ξ = τ/Ti ∈ [0, 1]. It is noted here that when p = 1 and if the shape functions for a linear beam element are used, then Eqs. (8) and (9) will give the element matrices of the TFE in [14]. 3. A simplified variant of time finite element: the ap-TFE In Eqs. (8) and (9), for an n-DOF system, the matrix H used for bp-TFEs is of the form    H1 H2 H6     H1 H2 H6   ···    H =  . . .  (11)    H1 H2 H6   | {z } | {z } | {z }  diag.matrix of size n diag.matrix of size n diag.matrix of size n The six shape functions here are exactly the ones for a compatible second-order finite beam with three equidistantly located nodes. With these shape functions, we can interpolate the displacement field and velocity field of a time element through the six nodal values of responses [15]. The size of H matrix is n × 6n, thus for bp-TFEs, the size of element “stiffness” matrix is 6n × 6n, and that of element “equivalent force” is 6n × 1. It raises to the idea of using non-compatible (spatial) finite axial bar elements from which the displacement field and velocity field can still be interpolated with polynomials of high degree through a lower number of nodal values, only at the beginning and at the end of the time interval (see Fig.1). For each TFE now, we select only the displacements and velocities at its two ends as the unknown nodal values of responses. Following this idea, the number of shape functions needed to interpolate the fields of responses of the time element can be reduced to 4, so that the matrix H now becomes of lower size n × 4n as    H1 H2 H4     H1 H2 H4   ···    H =  . . .  (12)    H1 H2 H4   | {z } | {z } | {z }  diag.matrix of size n diag.matrix of size n diag.matrix of size n We have plenty of choice for shape functions of non-compatible axial bar elements. Since the formulation shown in Eq. (8) requires second derivatives of H, in this article, the shape functions are chosen as polynomials of degree three as follows 2 3  2 3 H1 (ξ) = 1 − 3ξ + 2ξ , H2 (ξ) = ξ − 2ξ + ξ Ti (13) 2 3  2 3 H3 (ξ) = 3ξ − 2ξ , H4 (ξ) = −ξ + ξ Ti (14) 45
5. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering where ξ = τ/Ti ∈ [0, 1], and Ti is the “length” of the ith element. Derivation for the above shape functionsthethe above follows above shapeshape the procedure functions functions that followsfollows is commonly thethe procedure procedure used in the that that conventional is is commonly commonly FEM used used to obtain in in the the the shape functionsconventionalconventional for high-order FEMFEM to incompatibleto obtain obtain the the shape shape elements. functions functions With for thisfor high high choice-order-order of incompatible shapeincompatible functions, elements elements the size. . of el- ementWithWith “stiffness” thisthis choicechoice matrix of of shape shape can reduce functions, functions, to 4 nthe the× size 4sizen instead of of element element of 6 “stiffness”n “stiffness”× 6n as inmatrix matrix [15]. can Incan addition,reduce reduce to to with this choice of shape functions, the size of element “equivalent force” vector can reduce to 4n × 1 instead 44��××44�� insteadinstead of of 6 6��××66�� as as in in [15]. [15]. In In addition addition, ,with with this this choice choice of of shape shape functions, functions, of 6n × 1. thethe sizesize ofof elementelement “equivalent “equivalent force” force” vector vector can can to to 4 4��××11 instead instead of of 6 6��××1.1 . (a) (b) Figure 1. (a) Conventional (spatial) incompatible axial bar finite element; and (b) Corresponding ap-TFE FigureFigure 1.1. (a)(a) Conventional Conventional (spatial) (spatial) incompatible incompatible axial axial bar bar finite finite element; element; and and (b) (b) CorrespondingCorresponding ap ap-TFE-TFE When all element matrices are already determined, the global stiffness matrix and global equiv- alent force vector can be found by usual assembling procedure as in conventional finite element method,WhenWhen as all all element element matrices matrices are are already already determined, determined, the the global global stiffness stiffness matrix matrix and and global global Ksysqsys = fsys (15) equivalentequivalent forceforce vector vector can can be be found found by by usual usual assembling assembling procedure procedure as as in in conventional conventional sys wherefinitefiniteq elementcollectselement method allmethod nodal, ,as as responses in all temporal elements of the system. Eq. (15) can be effi- ciently solved by several usual algorithms, including the parallel method [14]. ��syssys��syssys==�sys�sys (15(15) ) sys 4. Numericalwherewhere ��sys examples collectscollects all all nodal nodal responses responses in in all all temporal temporal elements elements of of the the system. system. Eq. Eq. (15) (15) cancan bebe efficiently efficiently solved solved by by several several usual usual algorithms, algorithms, including including the the parallel parallel method method [1 [14].4 ]. 4.1. Forced vibration of an undamped SDOF system 44 NumericalNumerical examples examples Given a SDOF system having the mass of m = 1, the stiffness of k = π2/4, and there is no damping4.1.4.1. c ForcedForced= 0. The vibrationvibration system of of is an atan restundamped undamped when itSDOF SDOF is suddenly system system enforced by a pulse force f (t) given below GivenGiven aa SDOFSDOF systemsystem having having the the mass mass of of � �==11, ,the the stiffness stiffness of of � �==��/4/,4 and, and there there 1, 0 < t < 1 isis nono dampingdamping � �==00. .The The system system is fis at( tat) rest =rest when when it it is is suddenly suddenly enforced enforced by by a apulse pulse force force 0, t ≥ 1 ��((��)) givengiven belowbelow The exact displacement and velocity responses11, , 00 of<< the��<< system11 are found to be ��((��))== 00, , ��≥≥11  2 4 (1 − cos (πt/2)) /π , 0 < t < 1 TheThe exactexact displacementdisplacementu (t) and =and velocity velocity responses responses of of the the system system are are found found to to be be 4 (sin (πt/2) − cos (πt/2)) /π2, t ≥ 1 4(1 − cos(��/2))/�, 0 < � < 1 ( ) 4(1 − cos(��/2))/� , 0 < � < 1 ��(��)== and 44((sinsin((����//22))−−coscos((����//22))/)/��,, ��≥≥11  and and 2 sin (πt/2) /π, 0 < t < 1 u˙ (t) =  2sin2 ((cos��/(2π)t//2�),+ sin (πt/2)) /π, 0t <≥ �1 < 1 �̇ (�) = 2 sin(��/2) /�, 0 < � < 1 �̇ (�) = 2(cos(��/2) + sin(��/2))/�, � ≥ 1 For comparison, the problem2( iscos also(�� solved/2) + bysinβ(-Newmark��/2))/�, method� ≥ 1 - the most popular numerical integration method for transient structural dynamics - with three different sets of parameters (γ, β), 46
6. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering namely (0, 0), (1/2, 1/4), and (1/2, 1/6). The most accurate results obtained by using β-Newmark methods with these sets of parameters are shown in Tables1 and2, together with the results offered by [14, 15]. In [15], among results from bp-TFEs with different values of p, only the best result and the result with p = 1 are listed. The responses over the time interval [0, T] are considered, where T = 12 seconds. The time inter- val [0, T] is uniformly divided into N elements such that mesh size is τ = T/N. Four different mesh sizes are considered, namely τ = 1/2; 1/4; 1/8; and 1/16. For each mesh size, to quantify specifi- cally the accuracy, discrete max-norm errors, defined as maximum of absolute point-wise errors at all time nodes are computed and listed in Tables1 and2. The L2 norm error for displacement response 2 ku − uτk0 and the L norm error for velocity response ku˙ − u˙τk0 are plotted in Fig.2. It is seen that these L2 norm errors possess a quadratic convergence rate, a bit lower than that from [14]. The lower Table 1. Discrete max-norm error in displacement with the use of Newmark (γ, β) method and with the use of variants of TFEs Mesh size τ 1/2 1/4 1/8 1/16 Newmark (γ, β) 2.0725 · 10−1 5.6488 · 10−2 1.4457 · 10−2 3.6345 · 10−2 From [14] 5.6461 · 10−4 6.0392 · 10−5 1.0574 · 10−5 2.3831 · 10−6 bp-TFEs Best 7.5843 · 10−6 4.8081 · 10−7 2.9627 · 10−8 1.8455 · 10−9 [15] p = 1 3.7165 · 10−5 2.2662 · 10−6 1.4023 · 10−7 8.7608 · 10−9 p = −1 3.0051 · 10−3 7.7250 · 10−4 1.9172 · 10−4 4.7892 · 10−5 p = 0 4.2177 · 10−3 1.0449 · 10−3 2.6062 · 10−4 6.5117 · 10−5 ap-TFEs p = 1 1.5075 · 10−2 3.6228 · 10−3 8.9415 · 10−4 2.2342 · 10−4 p = 2 4.2525 · 10−2 9.4843 · 10−3 2.2943 · 10−3 5.6710 · 10−4 p = 3 1.2681 · 10−1 2.5448 · 10−2 5.9493 · 10−3 1.4609 · 10−3 Table 2. Discrete max-norm error in velocity with the use of Newmark (γ, β) method and with the use of variants of TFEs Mesh size τ 1/2 1/4 1/8 1/16 Newmark (γ, β) 3.4703 · 10−1 9.7632 · 10−2 2.5006 · 10−2 6.2638 · 10−3 From [14] 1.0041 · 10−2 2.4368 · 10−3 6.0898 · 10−4 1.5183 · 10−4 Best 3.1248 · 10−5 1.9402 · 10−6 1.2238 · 10−7 7.6559 · 10−9 bp-TFEs p = 1 1.0278 · 10−4 6.2780 · 10−6 3.9268 · 10−7 2.4459 · 10−8 p = −1 1.2759 · 10−2 3.1953 · 10−3 8.0830 · 10−4 2.0235 · 10−4 p = 0 2.0884 · 10−2 5.3131 · 10−3 1.3232 · 10−3 3.3096 · 10−4 ap-TFEs p = 1 4.2624 · 10−2 1.0111 · 10−2 2.5148 · 10−3 6.2557 · 10−4 p = 2 1.0224 · 10−1 2.2706 · 10−2 5.4982 · 10−3 1.3621 · 10−3 p = 3 2.8097 · 10−1 5.6330 · 10−2 1.3178 · 10−2 3.2294 · 10−3 47
7. � = 1 1.5075 ⋅ 10 3.6228 ⋅ 10 8.9415 ⋅ 10 2.2342 ⋅ 10 � = 2 4.2525 ⋅ 10 9.4843 ⋅ 10 2.2943 ⋅ 10 5.6710 ⋅ 10 � = 3 1.2681 ⋅ 10 2.5448 ⋅ 10 5.9493 ⋅ 10 1.4609 ⋅ 10 Table 2. Discrete max-norm error in velocity with the use of Newmark (�, �) method; and with the use of variants of TFEs Mesh size � 1/2 1/4 1/8 1/16 Newmark (�, �) 3.4703 ⋅ 10 9.7632 ⋅ 10 2.5006 ⋅ 10 6.2638 ⋅ 10 From [14] 1.0041 ⋅ 10 2.4368 ⋅ 10 6.0898 ⋅ 10 1.5183 ⋅ 10 bp- Best 3.1248 ⋅ 10 1.9402 ⋅ 10 1.2238 ⋅ 10 7.6559 ⋅ 10 TFEs � = 1 1.0278 ⋅ 10 6.2780 ⋅ 10 3.9268 ⋅ 10 2.4459 ⋅ 10 � = −1 1.2759 ⋅ 10 3.1953 ⋅ 10 8.0830 ⋅ 10 2.0235 ⋅ 10 � = 0 2.0884 ⋅ 10 5.3131 ⋅ 10 1.3232 ⋅ 10 3.3096 ⋅ 10 TFEs � = 1 4.2624 ⋅ 10 1.0111 ⋅ 10 2.5148 ⋅ 10 6.2557 ⋅ 10 - ap � = 2 1.0224 ⋅ 10 2.2706 ⋅ 10 5.4982 ⋅ 10 1.3621 ⋅ 10 Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering � = 3 2.8097 ⋅ 10 5.6330 ⋅ 10 1.3178 ⋅ 10 3.2294 ⋅ 10 convergence rate when compared to that from [14] is not always the case, but is case-dependent. Figure 2. Convergence for an undamped SDOF system As observed from Tables1 and2, the smaller the mesh size is, the higher magnitude order of accuracy. The ap-TFEs give comparable, though a little bit less, accurate results when compared even to the results offered by [14], not to mention about the comparison with bp-TFEs. However, when compared with the results offered by the best method among the β-Newmark family, the results from ap-TFEs are more accurate, even for the worst case. There is also an interesting observation. Similar to the results from bp-TFEs, although the cases with p = −1 or p = 0 are not supported by strictly mathematical arguments, the corresponding results are still more accurate than that obtained with ap-TFEs with higher values of p. 4.2. Forced vibration of a damped SDOF system The system in the first example is now changed with the damping property added. The damping coefficient considered is c = 0.2π. The initial conditions are u (0) = 4/π2 and u˙ (0) = 0. The exact responses of the system are found to be  4/π2, 0 < t < 1 u (t) =   −ζωnt e (A cos (ωDt) + B sin (ωDt)) , t ≥ 1 and  0, 0 < t < 1    u˙ (t) =  A (−ζωn cos (ωDt) − ωD sin (ωDt)) +  e−ζωnt   , t ≥    1  B (−ζωn sin (ωDt) + ωD cos (ωDt)) q p 2 where ωn = k/m, ζ = c/ (2mωn), ωD = ωn 1 − ζ , and A, B are solution of the following equation   " # " #  cos ωD sin ωD  A 4eζωn /π2  q q  =  2 2   −ζ cos ωD − 1 − ζ sin ωD −ζ sin ωD + 1 − ζ cos ωD  B 0 48
8. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering Table 3. Discrete max-norm error in displacement with the use of Newmark (γ, β) method and with the use of variants of TFEs Mesh size τ 1 1/2 1/4 1/8 Newmark (γ, β) 6.2597 · 10−2 1.6248 · 10−2 4.2875 · 10−3 1.0858 · 10−3 From [14] 3.2618 · 10−3 7.6250 · 10−4 1.8740 · 10−4 4.7333 · 10−5 p = −1 3.1360 · 10−3 7.2839 · 10−4 1.7896 · 10−4 4.4597 · 10−5 p = 0 3.9057 · 10−3 9.1779 · 10−4 2.2595 · 10−4 5.7094 · 10−5 ap-TFEs p = 1 4.8398 · 10−3 1.1319 · 10−3 2.8783 · 10−4 7.1898 · 10−5 p = 2 6.9592 · 10−3 1.5771 · 10−3 4.0280 · 10−4 9.9441 · 10−5 p = 3 1.0113 · 10−2 2.4408 · 10−3 5.9590 · 10−4 1.4906 · 10−4 Table 4. Discrete max-norm error in velocity with the use of Newmark (γ, β) method and with the use of variants of TFEs Mesh size τ 1 1/2 1/4 1/8 Newmark (γ, β) 9.3062 · 10−2 2.5437 · 10−2 6.7757 · 10−3 1.7218 · 10−3 From [14] 2.3461 · 10−2 5.7848 · 10−3 1.5065 · 10−3 3.7517 · 10−4 p = −1 2.1751 · 10−2 5.3271 · 10−3 1.3971 · 10−3 3.4834 · 10−4 p = 0 2.4464 · 10−2 6.1412 · 10−3 1.5824 · 10−3 3.9410 · 10−4 ap-TFEs p = 1 2.7502 · 10−2 7.1079 · 10−3 1.7921 · 10−3 4.4844 · 10−4 p = 2 3.1272 · 10−2 8.3309 · 10−3 2.0461 · 10−3 5.1651 · 10−4 p = 3 3.6942 · 10−2 1.0083 · 10−2 2.4470 · 10−3 6.1067 · 10−4 As in Example 1, we consider four different mesh sizes, including a coarser mesh, namely τ = 1; 1/2; 1/4; and 1/8. For this problem, among three options of (γ, β) of (0, 0) , (1/2, 1/4), and (1/2, 1/6) in the β-Newmark, the method with γ = 1/2 and β = 1/6 gives most accurate results as listed in the Tables3 and4 for comparison. Again here, better accuracy given by ap-TFEs can be pointed out, compared with ones offered by β-Newmark methods. For comparison with results from [14], the ap- TFE with p = −1 (shown in bold numbers) gives higher accuracy. However, this is not true with other values of p where the accuracy offered by ap-TFEs are often lower. Pointwise errors in results for the case when τ = 1/8 are depicted in Fig.3 where a quick decay can be observed as time evolves. Also, for some interval of time, the ap-TFE can give better accuracy. However, the overall performance in most of the time considered of ap-TFEs with p = 1 is worse than that from [14]. Fig.4 shows the convergence rate of the errors which is also approximately equal to 2. 49
9. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering Figure 3. History of pointwise errors of responses under mesh � = 1/8 for a damped SDOF system Figure 3. HistoryFigure 3 of. History pointwise of pointwise errors of errors responses of responses under under mesh meshτ = 1 �/8=for1/8 a dampedfor a damped SDOF system SDOF system FigureFigure 4. 4. Convergence Convergence forfor aa damped SDOF SDOF system system Figure 4. Convergence for a damped SDOF system 4.3. Free vibration of damped MDOF system A 2-DOF system subjected to an initial condition is shown in Fig.5. The total mass of the bar is m = 1. The spring stiffness is k = 1, and the damping coefficient of the dampers is c = 0.2. With the above values, the structural matrices can be obtained as " # " # " # 1 2 1 1 0 1 0 M = × , C = 0.2 × , K = 6 1 2 0 1 0 1 h iT h iT The initial conditions are u (0) = 1.0 0.0 and u˙ (0) = 0.0 0.0 . The exact responses of the system are found to be [14] 50
10. 4.3. Free vibration of damped MDOF system A 2-DOF system subjected to an initial condition is shown in Figure 5. The total mass of the bar is � = 1. The spring stiffness is � = 1, and the damping coefficient of the dampers is � = 0.2. With the above values, the structural matrices can be obtained as Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering 1 2 1 1 0 1 0 � = × , � = 0.2 × , � = 6 1 2 0 1 0 1 The initial conditions are �(0) = [1.0 0.0] and �̇ (0) = [0.0 0.0]. The exact 2   1 X responses of the system areζ found to be [14] −ζiωit  i  u (t) = φie cos ωDit + q sin ωDit 2  2  i=1  1 − ζ  � i !! �(�) = �� cos �� + sin �� 2 1 − � X2 ω 1 i −ζiωit u˙ (t) = − φi q e sin ωDit 2 2 1 � i=1 − !! 1 ζi �̇ (�) = − � � sin �� 2 1 − � p p h iT h iT with ω1 = 6k/m, ω2 = 2k/m, ζ1 = 3withc/ (m �ω 1=),ζ62�=/�c/, �(m ω=2),2φ�1/�=, � 1= 3−�/1(��,φ),2 �== �1/(� 1�), ,� = [1 −1] , q 2 � = [1 1] , and � = �1 − � . and ωDi = ωi 1 − ζi . Uniform meshes are adopted with mesh sizes τ = 1/2, 1/4, 1/8, and 1/16. The pointwise errors in the results obtained with these meshes using ap- TFEs of different values of p are shown in Tables5 and6, side-by-side with results from [14] and from the β-Newmark method with γ = 1/2, β = 1/6 in the finest mesh size τ = 1/16. It is noted here that, in [14], there is no data for the case when the FigureFigure 5. Free 5. Free vibration vibration of a ofdamped a damped 2-DOF 2-DOF system mesh size is τ = 1/8. Since the accuracyUniform given meshes are adopted with meshsystem sizes � = 1/2, 1/4,1/8, and 1/16. The by bp-TFEs is much higher than all otherpointwise variants errors in the results obtained with these meshes using ap-TFEs of different of TFEs, as shown in Example 1, we dovalues not of include � are shown theresults in Table from 5 andusing Table 6, bp-TFEs side-by-side here with for results the from [14] and comparison. In these tables, the numbersfrom in the bold �-Newmark format method show thewith cases� = 1/ in2, � which= 1/6 the in the use finest of mesh ap- size � = 1/16. It is noted here that, in [14], there is no data for the case when the mesh size is � = 1/8. Since the accuracy given by bp-TFEs is much higher than all other variants of TFEs, as Table 5. Discrete max-norm error in displacements with the use of Newmark (γ, β) method shown in Example 1, we do not include the results from using bp-TFEs here for the and with the use of variants of TFEs comparison. In these tables, the numbers in bold format show the cases in which the use of ap-TFEs gives better accuracy when compared with the corresponding values from Mesh size τ 1/2 1/4 1/8 1/16 [14]. This finding is supportive for the comment we made following Example 1 saying Newmark (γ, β) that the conventional TFEs do not always give better6 accuracy.3087 · 10against−4 the ap-TFEs. The From [14] 3.9984 · 10−3 9.5821 · 10−4 - 6.0209 · 10−5 p = −1 3.9079 · 10−3 9.8772 · 10−4 2.5021 · 10−4 6.2460 · 10−5 u1 p = 0 4.3529 · 10−3 1.0445 · 10−3 2.6396 · 10−4 6.5828 · 10−5 ap-TFEs p = 1 4.8662 · 10−3 1.1483 · 10−3 2.8287 · 10−4 7.1029 · 10−5 p = 2 5.5823 · 10−3 1.2968 · 10−3 3.1816 · 10−4 7.9125 · 10−5 p = 3 6.8782 · 10−3 1.5723 · 10−3 3.8374 · 10−4 9.5234 · 10−5 Newmark (γ, β) 9.4162 · 10−4 From [14] 3.0197 · 10−3 7.1961 · 10−4 - 4.4264 · 10−5 p = −1 3.0519 · 10−3 7.3939 · 10−4 1.8471 · 10−4 4.6400 · 10−5 u2 p = 0 3.0428 · 10−3 7.2326 · 10−4 1.7840 · 10−4 4.4449 · 10−5 ap-TFEs p = 1 2.9384 · 10−3 6.8345 · 10−4 1.6789 · 10−4 4.1912 · 10−5 p = 2 2.7087 · 10−3 6.8871 · 10−4 1.6666 · 10−4 4.1229 · 10−5 p = 3 4.7642 · 10−3 1.1544 · 10−3 2.7986 · 10−4 6.9635 · 10−5 51
11. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering Table 6. Discrete max-norm error in velocities with the use of Newmark (γ, β) method and with the use of variants of TFEs Mesh size τ 1/2 1/4 1/8 1/16 Newmark (γ, β) 1.2609 · 10−3 From [14] 2.5823 · 10−2 6.2776 · 10−3 - 3.9686 · 10−4 p = −1 2.4327 · 10−2 5.9301 · 10−3 1.4988 · 10−3 3.7484 · 10−4 u˙1 p = 0 2.6025 · 10−2 6.3236 · 10−3 1.6020 · 10−3 3.9987 · 10−4 ap-TFEs p = 1 2.7866 · 10−2 6.7484 · 10−3 1.7130 · 10−3 4.2720 · 10−4 p = 2 3.0016 · 10−2 7.2427 · 10−3 1.8435 · 10−3 4.5952 · 10−4 p = 3 3.2989 · 10−2 7.9196 · 10−3 2.0280 · 10−3 5.0509 · 10−4 Newmark (γ, β) 1.9733 · 10−3 From [14] 1.8854 · 10−2 4.5497 · 10−3 - 2.8106 · 10−4 p = −1 1.7774 · 10−2 4.3001 · 10−3 1.0662 · 10−3 2.6603 · 10−4 u˙2 p = 0 1.8856 · 10−2 4.5463 · 10−3 1.1259 · 10−3 2.8080 · 10−4 ap-TFEs p = 1 1.9832 · 10−2 4.7654 · 10−3 1.1789 · 10−3 2.9390 · 10−4 p = 2 2.0503 · 10−2 4.9123 · 10−3 1.2143 · 10−3 3.0267 · 10−4 p = 3 2.1167 · 10−2 5.7390 · 10−3 1.4071 · 10−3 3.4948 · 10−4 TFEs gives better accuracy when compared with the corresponding values from [14]. This finding is supportive for the comment we made following Example 1 saying that the conventional TFEs do not always give better accuracy against the ap-TFEs. The responses computed with the mesh size τ = 1/2 responses computed with the mesh size � = 1/2 are depicted in Figure 6 for schematic are depicted inview. Fig.6 for schematic view. FigureFigure 6. Responses6. Responses for for aa dampeddamped rigid rigid bar bar under under mesh mesh � =τ1=/21 /2 Table 5. Discrete max-norm error in displacements with the use of Newmark (�, �) method; and with the use of variants of TFEs Mesh size � 1/2 521/4 1/8 1/16 Newmark (�, �) - - - 6.3087 ⋅ 10 From [14] 3.9984 ⋅ 10 9.5821 ⋅ 10 - 6.0209 ⋅ 10 � = −1 �. ���� ⋅ ��� 9.8772 ⋅ 10 2.5021 ⋅ 10 6.2460 ⋅ 10 � � = 0 4.3529 ⋅ 10 1.0445 ⋅ 10 2.6396 ⋅ 10 6.5828 ⋅ 10 TFEs � = 1 4.8662 ⋅ 10 1.1483 ⋅ 10 2.8287 ⋅ 10 7.1029 ⋅ 10 - ap � = 2 5.5823 ⋅ 10 1.2968 ⋅ 10 3.1816 ⋅ 10 7.9125 ⋅ 10 � = 3 6.8782 ⋅ 10 1.5723 ⋅ 10 3.8374 ⋅ 10 9.5234 ⋅ 10 Newmark (�, �) - - - 9.4162 ⋅ 10 � From [14] 3.0197 ⋅ 10 7.1961 ⋅ 10 - 4.4264 ⋅ 10 - s ap � = −1 3.0519 ⋅ 10 7.3939 ⋅ 10 1.8471 ⋅ 10 4.6400 ⋅ 10 TFE
12. Nguyen, T. X., Tran, L. T. / Journal of Science and Technology in Civil Engineering 5. Conclusions Motivated by effective use of high-order time finite elements in solving dynamic problems shown in [14, 15], a modification in choosing shape functions was made to have a simplified variant of time finite elements, denoted by ap-TFE. In this article, the ap-TFEs are obtained using four shape functions, opposed to six ones as for bp-TFEs, in the form of polynomials of degree three, opposed to degree five as for bp-TFEs. With this proposal, the size of matrix H is reduced to n × 4n, instead of n × 6n as in the case of bp-TFEs, leading to the size of element “stiffness” matrix of 4n × 4n, instead of 6n × 6n as in the case of bp-TFEs. As a compromise, the accuracy given by using ap-TFEs is much lower when compared with ones given by using bp-TFEs as shown in [15]. Nevertheless, the proposed ap-TFEs still give much better accuracy when compared with methods in β-Newmark family. They also possess a comparable accuracy and convergence rate when compared with the ones provided by [14]. The good performance of ap-TFEs illustrated through examples in this article sounds that this type of TFE would be an alternative effective tool for structural dynamic analysis. References [1] Zienkiewicz, O. C. (1977). A new look at the Newmark, Houbolt and other time stepping formulas. A weighted residual approach. Earthquake Engineering & Structural Dynamics, 5(4):413–418. [2] Newmark, N. M. (1962). A Method of Computation for Structural Dynamics. Transactions of the Ameri- can Society of Civil Engineers, 127(1):1406–1433. [3] Dokainish, M., Subbaraj, K. (1989). A survey of direct time-integration methods in computational struc- tural dynamics—I. Explicit methods. Computers & Structures, 32(6):1371–1386. [4] Subbaraj, K., Dokainish, M. A. (1989). A survey of direct time-integration methods in computational structural dynamics—II. Implicit methods. Computers & Structures, 32(6):1387–1401. [5] Hughes, T. J. R., Belytschko, T. (1983). A Pre´cis of Developments in Computational Methods for Tran- sient Analysis. Journal of Applied Mechanics, 50(4b):1033–1041. [6] Park, S. (1997). Development and applications of finite elements in time domain. Doctoral dissertation, Virginia Polytechnic Institute and State University. [7] Argyris, J. H., Scharpf, D. W. (1969). Finite Elements in Time and Space. The Aeronautical Journal, 73 (708):1041–1044. [8] Simkins, T. E. (1978). Unconstrained Variational Statements for Initial and Boundary-Value Problems. AIAA Journal, 16(6):559–563. [9] Fried, I. (1969). Finite-element analysis of time-dependent phenomena. AIAA Journal, 7(6):1170–1173. [10] Zienkiewicz, O. C., Parekh, C. J. (1970). Transient field problems: Two-dimensional and three- dimensional analysis by isoparametric finite elements. International Journal for Numerical Methods in Engineering, 2(1):61–71. [11] Hulbert, G. M. (1992). Time finite element methods for structural dynamics. International Journal for Numerical Methods in Engineering, 33(2):307–331. [12] French, D. A., Peterson, T. E. (1996). A continuous space-time finite element method for the wave equation. Mathematics of Computation, 65(214):491–506. [13] Tang, W., Sun, Y. (2012). Time finite element methods: A unified framework for numerical discretizations of ODEs. Applied Mathematics and Computation, 219(4):2158–2179. [14] Wang, L., Zhong, H. (2017). A time finite element method for structural dynamics. Applied Mathematical Modelling, 41:445–461. [15] Nguyen, T. X., Tran, L. T. (2022). A High-Order Time Finite Element Method Applied to Structural Dynamics Problems. In Modern Mechanics and Applications, Springer, 137–148. 53