- Email: [email protected]
Numerical methods solving finite element dynamic equations
Appendix
Numerical methods solving finite element dynamic equations This appendix gives some numerical methods, developed by Xing (2005), in FORTRAN designed to solve dynamic finite element (FE) equations. The principles, solution processes, examples, and comparison of five time integration methods—Newmark, Wilson-θ, Hilber-α, Hilber collocation, and central difference—which are often used to solve dynamic FE equations, are presented. A computer code Program of Time Element Methods (PTEM), developed by the author, is introduced and presented. This code provides a generalized functional module using these five algorithms to solve the FE equation t M x€ 1 t C x_ 1 t Kx 5 Rðt; x; x_ Þ describing linear and nonlinear structural dynamics. This functional program module is conveniently incorporated into any computer program used to investigate complex dynamics problems in maritime engineering.
A.1
Introduction
Direct integration methods, or time element methods, are important and necessary approaches in analyzing dynamic response problems of linear and nonlinear dynamic systems (Bathe and Wilson, 1976; Bathe, 1982, 1986, 1996; Belytschko and Hughes 1983; Zheng and Tan, 1985; Zheng 1986). Essentially, direct numerical integration methods are based on two ideas. First, instead of satisfying the FE dynamic equation at any time t, these methods are aimed to satisfy it only at discrete time intervals Δt apart. The second is that a variation of displacements, velocities, and accelerations within each time interval is assumed. In depending on the different discrete points and variation assumptions used, there are different time integration schemes with different performances of accuracy, stability, and cost of solution. The central difference method is an often used explicit method with three-point interpolation. Houbolt (1950) developed a four-point implicit scheme. Newmark (1959) introduced two parameters γ and β to represent displacement and velocity, respectively, and created a very useful algorithm called the Newmark method. Hilber and Taylor (1977) introduced another parameter, α, in order to control numerical damping and to improve accuracy, into the Newmark scheme and created the Hilber-α method. Wilson et al. (1973), assuming a linear acceleration distribution over a time interval (t,t 1 θΔt), constructed the Wilson-θ method. Combining the ideas of the Wilson and Newmark methods, Hilber and Hughes (1978) developed another method called the Hilber collocation method. Zienkiewiez (1977) investigated the constructions of the Newmark and other time integration methods using a weighted residual approach. Based on variational principles and interpolations in the time domain, Xing and Zheng (1985, 1987) studied existing methods and created a new approach. In an integrated variational frame for the interpolations in space and time domains, the concept of a four-dimensional FE method was proposed. In the current FE software, the five methods—the central difference, Newmark, Wilson-θ, Hilber-α, Hilber collocation schemes—are widely used. This appendix gives a description and comparison of these five methods. A computer code PTEM, developed by the author, is introduced and listed here for application by readers.
A.2 A.2.1
Fundamental iteration formulations Generalized dynamic equation
In a linear FE analysis (FEA), a generalized dynamic equation is expressed as M€x 1 C_x 1 Kx 5 R;
(A.1)
577
578
Appendix: Numerical methods solving finite element dynamic equations
where M, C, and K are constant mass, damping, and stiffness matrices of the system, x, x_ , and x€ , represent the displacement, velocity, and acceleration vectors, and R is a time-dependent load vector. For nonlinear systems, Eq. (A.1) can be replaced by the following increment form: t
MΔ€x 1 t CΔ_x 1 t KΔx 5 Δt R;
(A.2)
where tM, tC, and tK denote the tangent mass, damping, and stiffness matrices at time t, and Δt R is the load increment. Since at each time step, the tangent mass, damping, and stiffness matrices are recalculated, the cost of solving Eq. (A.2) is quite high. For practical problems, nonlinearity is mainly caused by stiffness and damping but not involving the mass of the system. To reduce the computational cost, Eq. (A.2) is represented as t 1 Δt
M x€ ðkÞ 1 0 CΔ_xðkÞ 1 0 KΔxðkÞ 5 t 1 Δt R 2 t 1 Δt Fd ð_xðk21Þ Þ 2 t 1 Δt Fs ðxðk21Þ Þ;
t 1 Δt ðkÞ
5 t 1 Δt xðk 2 1Þ 1 ΔxðkÞ ;
t 1 Δt ðkÞ
5 t 1 Δt _x ðk 2 1Þ 1 Δ_xðkÞ ;
x
x_
(A.3)
k 5 1; 2; 3; :::
where 0C and 0K are the initial damping and stiffness matrices at time t 5 0, and Fd and Fs represent the damping and stiffness forces that can be expressed in a summation of linear and nonlinear components as follows: d Fd 5 0 C x_ 1 F~ ð_xÞ;
s Fs 5 0 Kx 1 F~ ðxÞ:
(A.4)
Eq. (A.3) is now rewritten in the form t 1 Δt
t 1 Δt ~ R 2 t 1 Δt F~ d ð_xðk21Þ Þ 2 t 1 Δt F~ s ðxðk21Þ Þ 5 t 1 Δt R; ðM€xðkÞ 1 C_xðkÞ 1 KxðkÞ Þ 5
(A.5)
where for convenience the superscripts of the damping and stiffness matrices are neglected. Eq. (A.5) provides an iteration formulation to obtain a solution x(k) at iteration k from previous result x(k21). The following three convergence criteria may be chosen to check whether the convergence is reached at iteration k.
A.2.2
Iteration criterions
A.2.2.1 Energy criterion ΔxðkÞ t 1 Δt ½R 2 Fd ð_xðk21Þ 2 Fs ðxðk21Þ 2 M€xðk21Þ Δxð1Þ U½t 1 Δt R 2 t Fd 2 t Fs 2 Mt1Δt x€ ð0Þ
(A.6)
# ETOL:
A.2.2.2 Force criterion :t1Δt ½R2Fd ð_xðk21Þ 2Fs ðxðk21Þ 2M€xðk21Þ :2 :½t1Δt R2 t Fd 2 t Fs 2Mt1Δt x€ ð0Þ :2
(A.7)
# RTOL:
A.2.2.3 Displacement criterion :ΔxðkÞ :2 # DTOL: :xðkÞ :2
Here :a:2 5
pffiffiffiffiffiffiffiffi aUa, ETOL, RTOL, and DTOL represent accuracy errors defined by users.
(A.8)
Appendix: Numerical methods solving finite element dynamic equations
A.3 A.3.1
579
Five numerical integration schemes Newmark method
The Newmark method assumes that the velocity and displacement vectors at time step n 1 1 take the forms x_ n11 5 x_ n 1 ½ð1 2 γÞ€xn 1 γ x€ n11 ; xn11 5 xn 1 x_ n Δt 1 ½ð0:5 2 βÞ€xn 1 β x€ n11 Δt2 ;
(A.9)
~ n11 ; M€xn11 1 C_xn11 1 Kxn11 5 R from which the displacement, velocity, and acceleration vectors xn11, x_ n11 , and x€ n11 at time step n 1 1 are derived from their values at the previous time step n. The unconditional stable condition of the Newmark method requires that β $ 0:25ð0:51γÞ2 :
γ $ 0:5;
A.3.2
(A.10)
Wilson-θ method
This method assumes that the acceleration in time interval τAðt; t 1 θΔtÞ varies in a linear form: τ ð€xn1θΔt 2 x€ n Þ: x€ n1τ 5 x€ n 1 θΔt
(A.11)
The integration of this equation gives that x_ n1τ 5 x_ n 1 x€ n τ 1
τ2 ð€xn1θΔt 2 x€ n Þ; 2θΔt
τ2 τ3 xn1τ 5 xn 1 x_ n τ 1 x€ n 1 ð€xn1θΔt 2 x€ n Þ; 2 6θΔt
(A.12)
which, at τ 5 θΔt, x_ n1θΔt , and x€ n1θΔt , can be solved in terms of xn1θΔt . The substitution of the results of x_ n1θΔt and x€ n1θΔt into the dynamic equation M€xn1θΔt 1 C_xn1θΔt 1 Kxn1θΔt ~ n1θΔt 5 R ~ n 1 θðR ~ n1Δt 2 R ~ nÞ 5R
(A.13)
provides a solution for xn1θΔt . The solution of xn1Δt , x_ n1Δt , and x€ n1Δt can be obtained by 3 6 6 x€ n 2 2 x_ n 1 3 2 ðxn1θΔt 2 xn Þ; x€ n11 5 1 2 θ θ Δt θ Δt x_ n11 5 x_ n 1
(A.14a)
Δt ð€xn11 1 x€ n Þ; 2
Δt2 ð€xn11 1 2€xn Þ: xn11 5 xn 1 Δtx_ n 1 6
(A.14b)
The unconditional stable condition of the method requires θ $ 1:37.
A.3.3
Hilber-α method
Based on the Newmark method, Hilber introduced a parameter α into Eq. (A.9) and obtained that ~ n11 ; M€xn11 1 C_xn11 1 K½xn11 ð1 1 αÞ 2 αxn 5 R
(A.15)
580
Appendix: Numerical methods solving finite element dynamic equations
developing a modification algorithm. The parameter α controls human-made damping and increases the accuracy. The unconditional stable condition of this method requires that γ 5 0:5 2 α; β 5 0:25ð1 2 αÞ;
A.3.4
20:5 # α # 0 :
(A.16)
Hilber collocation method
Combining the Wilson-θ and Newmark methods, Hilber constructed a collocation scheme using the following form: x€ n1θ 5 x€ n 1 θð€xn11 2 x€ n Þ; x_ n1θ 5 x_ n 1 θΔt½ð1 2 γÞ€xn 1 γ x€ n1θ ; xn1θ 5 xn 1 θΔtx_ n 1 ðθΔtÞ2 ½ð0:5 2 βÞ€xn 1 β x€ n1θ ; ~ n1θ ; M€xn1θ 5 C€xn1θ 1 Kxn1θ 5 R
(A.17)
in which x_ n1θ and x€ n1θ are represented in terms of xn1θ , and then xn1θ is obtained by solving ^ ^ n1θ 5 R; Kx 1 γ ^5 C 1 K; M1 K βθΔt βðθΔtÞ2 ^ 5 ð1 2 θÞR ~ n 1 θR ~ n11 ; R 0 1 3 2 1 1 1 xn 1 1 M4 x_ n 1 @ 2 1Ax€ n 5; βðθΔtÞ 2β βðθΔtÞ2 0 1 0 1 2 3 γ γ γ xn 1 @ 2 1Ax_ n 1 @ 2 1AθΔtx€ n 5: 1 C4 βðθΔtÞ β 2β
(A.18)
Hilber demonstrated that this method is unconditionally stable provided that γ 5 0;
θ 2θ2 2 1 $β$ ; 2ðθ 1 1Þ 4ð2θ3 2 1Þ
θ$1
(A.19)
This method reduces to a Newmark method of γ 5 0:5 and β 5 0:25 if θ 5 1 and to a Wilson-θ method if γ 5 1=2 and β 5 1=6.
A.3.5
Central difference method
For the central difference method, it is assumed that
x_ n 5
1 ðxn11 2 xn21 Þ; 2Δt (A.20a)
1 x€ n 5 2 ðxn11 2 2xn 1 xn21 Þ; Δt ~ n; M€xn 1 C_xn 1 Kxn 5 R
(A.20b)
Appendix: Numerical methods solving finite element dynamic equations
581
from which it follows that ^ n11 5 R; ^ Kx ^ 5 1 M 1 1 C; K Δt2 2Δt 0 1 0 1 2 1 1 ^ 5R ~ n 2 @K 2 CAxn21 ; R M A xn 2 @ 2 M 2 Δt2 Δt 2Δt x21 5 x0 2 Δtx_ 0 1
(A.21)
Δt2 x€ 0 : 2
The central difference method is an explicit method for which iteration is not needed, and therefore it is widely used in nonlinear analysis. This method is conditionally stable, and the stable condition requires Δt #
Tm ; π
(A.22)
where Tm is the minimum vibration period of the problem.
A.4
Stability and accuracy
As previously described, the four methods, except for the central difference method, are unconditionally stable, implying that for any time step Δt, an approximate solution can be definitely obtained. However, for a large time step Δt, the accuracy of the result is lower but at less of a cost of calculation. Here, the stability and accuracy analysis of the methods are given for a reference in applications.
A.4.1
Stability and spectral radius
The stability of an integration method is determined by examining the behavior of the numerical solution for arbitrary initial conditions. Assume that the iteration formulation of a method to determine the solution t 1 Δt x^ from the solution t x^ is defined as t 1 Δt
x^ 5 At x^ ;
(A.23)
where A represents the matrix operator of a method. Let ρ(A) be the spectral radius of A defined as ρðAÞ 5 max jλj j;
j 5 1; 2; 3; . . . :
(A.24)
The necessary and sufficient condition of a stable time integration method is ρðAÞ # 1:
A.4.2
(A.25)
Period elongation and amplitude decay
Various studies of the accuracies of time integration methods have been reported. Here, we only summarize some important solution characteristics. For this purpose, we consider the solution of the initial value problem defined by x€ 1 ω2 x 5 0;
0
x 5 1:0; 0 x_ 5 0:0;
(A.26)
for which the exact solution is x 5 cos ωt. However, there is a period elongation PE and an amplitude decay AD in the approximate solution obtained by time integration methods. It has been defined that Percentage amplitude decay 5 AD 3 100%; Percentage period elongation 5
PE 3 100%; T
(A.27)
582
Appendix: Numerical methods solving finite element dynamic equations
to describe the accuracy of a method. Also, a damping factor η has been introduced as AD 5 1 2 e22πη ;
η5
2lnð1 2 ADÞ : 2π
(A.28)
For the often used time integration methods, the references by Hilber et al. (1977) and Hilber and Hughes (1978) present the spectral radii ρ(A), the percentage amplitude decay, the damping factor, and the percentage of period elongation as functions of Δt/T. These curves provide a good comparison of the discussed time integration methods. From the point of view of engineering applications, all of the five methods can provide sufficient accuracy for engineering analysis, provided a suitable time step chosen.
A.5
Implementation
The implementation process of these five time integration methods is described as follows.
A.5.1
Initial calculations
1. Form mass, stiffness, and damping matrices M, K, and C. 2. Initialize 0 x, 0 x_ , and 0 x€ 5 M21 ð0 R 2 C0 x_ 2 K0 xÞ. For the central difference method, you need to calculate 21 x 5 0 x 5 0 x_ Δt 1 0:50 x€ Δt2 . 3. Select time step Δt , Δtcr , and calculate the integration constants.
A.5.1.1 Newmark method γ $ 0:5; Ak 5 1;
β $ 0:25ð0:51γÞ2 ; A0 5 1ðβΔt2 Þ; A1 5 γ=ðβΔtÞ;
A2 5 1=ðβΔtÞ;
A3 5 1ð2βÞ 2 1;
A5 5 ðγ=β 2 2ÞΔt=2;
A4 5 γ=β 2 1;
A6 5 ð1 2 γÞΔt;
A7 5 γΔt:
A.5.1.2 Hilber-α method 0:5 # α # 0;
γ 5 0:5 2 α; β 5 0:25ð12αÞ2 ;
Ak 5 1 1 α;
A0 5 1=ðβΔt2 Þ;
A2 5 1=ðβΔtÞ;
A1 5 γ=ðβΔtÞ;
A3 5 1=ð2βÞ 2 1;
A5 5 ðγ=β 2 2ÞΔt=2;
A4 5 γ=β 2 1;
A6 5 ð1 2 γÞΔt;
A7 5 γΔt:
A.5.1.3 Wilson-θ method θ $ 1:37;
Ak 5 1;
A0 5 6=ðθΔt2 Þ; A3 5 θΔt=2; A6 5 1 2 3=θ;
A1 5 3=ðθΔtÞ; A2 5 2A1 ;
A4 5 A0 =θ;
A5 5 2 A2 =θ;
A7 5 Δt=2; A8 5 Δt2 =6:
Appendix: Numerical methods solving finite element dynamic equations
583
A.5.1.4 Hilber collocation method θ $ 1;
γ 5 0:5;
Ak 5 1;
θ 2θ2 2 1 $β$ ; 2ðθ 1 1Þ 4ð2θ3 2 1Þ
A0 5 1=ðβθ2 Δt2 Þ;
A1 5 γ=ðβθΔtÞ;
A2 5 1=ðβθΔtÞ; A3 5 1=ð2βÞ 2 1; A4 5 γ=β 2 1; γ A0 21 2 1 θΔt; A6 5 A5 5 ; A7 5 2 ; 2β θ βθ Δt A8 5 2
1 ; 2βθ
B1 5 Δtγ;
A9 5 ð1 2 γÞΔt;
B2 5 ð0:5 2 βÞθ2 Δt2 ;
B3 5 βθ2 Δt2 :
A.5.1.5 Central difference method Δt # Tn =π 5 Δtcr ; Ak 5 0; A0 5 1=Δt2 ;
A.5.2
A1 5 1ð2ΔtÞ;
A2 5 2A0 :
Form effective stiffness matrix ^ 5 A0 M 1 A1 C 1 Ak K: K
A.5.3
Triangular decomposition of matrix K^ ^ 5 LDLT : K
A.5.4
(A.29a)
(A.29b)
Calculations at each time step
1. Form effective loads at time t 1 Δt. For a convenience in equation expression, we define that t
R 5 MðA0 t x 1 A2 t x_ 1 A3 t x€ Þ 1 CðA1 t x 1 A4 t x_ 1 A5 t x€ Þ:
(A.30)
Newmark method t 1 Δt
^ 5 t 1 Δt R 1 t R: R
(A.31)
t 1 Δt
^ 5 t 1 Δt R 1 t R: R
(A.32)
Hilber-α method
Wilson-θ method t 1 Δt
^ 5 A2 t R 1 θt 1 Δt R 2 A2 Mt x€ R 2 Kðt x 1 A3 t x_ 1 A5 t x€ Þ 2 Cðt x 1 A4 t x€ Þ:
(A.33)
Hilber collocation method t 1 Δt
^ 5 ð1 2 θÞt R 1 θt 1 Δt R 1 t R: R
(A.34)
584
Appendix: Numerical methods solving finite element dynamic equations
Central difference method t 1 Δt
^ 5 t R 1 MðA2 t x 2 A0 t 2 Δt xÞ 2 Kt x 1 CA1 t 1 Δt x: R
(A.35)
2. Solve for x^ at time t 1 Δt: ^ ^ x 5 LDLT x^ 5 t 1 Δt R: K^
(A.36)
3. Evaluate displacements, velocities, and accelerations at time t 1 Δt. Newmark and Hilber-α methods t 1 Δt
x 5 x^ ;
t 1 Δt
x€ 5 A0 ð^x 2 t x^ Þ 2 A2 t x_ 2 A3 t x€ ;
(A.37)
x 5 t x_ 1 A6 t x€ 1 A7 t 1 Δt x€ :
t 1 Δt _
Wilson-θ method t 1 Δt
x€ 5 A4 ð^x 2 t x^ Þ 1 A5 t x_ 1 A6 t x€ ; x 5 t x_ 1 A7 ðt x€ 1 t 1 Δt x€ Þ;
t 1 Δt _ t 1 Δt
(A.38)
x 5 t x 1 Δtt x_ 1 A8 ð2t x€ 1 t 1 Δt x€ Þ:
Hilber collocation method t 1 Δt
x€ 5 A6 ð^x 2 t x^ Þ 1 A7 t x_ 1 A8 t x€ ; x 5 t x_ 1 A9 t x€ 1 B1 t 1 Δt x€ ;
t 1 Δt _ t 1 Δt
x 5 x 1 Δt x 1 B2 x€ 1 B3 t
t_
t
(A.39) t 1 Δt
x€ :
Central difference method t 1 Δt
x 5 x^ ;
x 5 A1 ðt 1 Δt x 2 t 2 Δt xÞ;
t_ t
A.6
(A.40)
x€ 5 A0 ðt 1 Δt x 2 2t x 1 t 2 Δt xÞ:
Time element program
The time element program PTEM, developed by Xing (1989), is a functional module that includes the five time integration methods just discussed. This program is used to solve linear and nonlinear dynamic FEA equations to obtain displacements, velocities, and accelerations of the system at time t 1 Δt from their values at time t. It can be incorporated into an FEA program as a functional module or directly used to solve dynamic equations obtained using other software.
A.6.1
Main characteristics
1. The program consists of about 1200 FORTRAN sentences. All functions of the program are completed in the module. It provides a convenient interface for users to connect the module into other FEA programs. 2. There are three forms to solve Eq. (A.1), which is controlled by a control variable JC as follows: JC 5 1: The matrices M, C, and K are assumed unchanged at each time step, and load R is dependent on time only. This function is mainly used to solve linear problems with no iteration required at each time step. JC 5 2: The matrices M, C, and K are still assumed unchanged at each time step, but the load R is dependent on unknown displacement and velocity. For the four time integration methods except the central difference method, iteration is required at each time step. Users can choose the convergence criteria. This process is mainly used to solve problems with weak nonlinearity. For problems with strong nonlinearity, iteration convergence may not reach, and the central difference method can be used.
Appendix: Numerical methods solving finite element dynamic equations
585
JC 5 3: The matrices M, C, and K are changed at each time step, but the load R is independent of only the time function unknown variables. During the calculation at each time step the matrices M, C, and K are required to be updated using FEA software. There are three flags JCM, JCC, and JCK to control the changes of matrices M, C, and K, respectively. 3. A control flag IJK is designed for users to choose one of the five time integration methods. 4. The program can automatically check the parameters α, β, γ, and θ to satisfy the stability conditions. 5. The load functions designed in the program include sinusoidal functions, discrete time functions given by input data, and the nonlinear spring function FUNC1 5 Ax2 1 Bx3
(A.41)
FUNC1 5 A_x2 1 B_x3 ;
(A.42)
and nonlinear damping function
where A and B are constants defined by users. Users can input their load functions, defined by themselves, and define the starting time of each load. 6. Output results are the displacements, velocities, and accelerations at the selected degrees of freedom (DOFs) by users. Users can choose the time intervals of the output data.
A.6.2
Flowchart of the program
A flowchart of the program is shown in Fig. A.1. All the input data, the values of control variables, and the method parameters are supplied by an INPUT.DAT file prepared by users. A file OUTPUT.DAT produced automatically by the program gives the output results.
A.7 A.7.1
Examples Example 1
A nondamping free vibration of a system with 1 DOF is defined by the following equation: x€ 1 x 5 0;
xð0Þ 5 1:0;
_ 5 1:0 ; xð0Þ
(A.43)
which has a theoretical solution x5
pffiffiffi π 2 cos t 5 : 4
(A.44)
Table A.1 lists the displacements calculated using the program in association with the parameters Δt 5 π=8, γ 5 0:5, β 5 0:25, α 5 20:01, and θ 5 1:4 for the Wilson-θ method and θ 5 1:001 for the Hilber collocation method.
A.7.2
Example 2
Let us consider a system with 2 DOFs described by the following equation: 2 4
2 0 0 1
2
54
x€1 x€2
3
2
514
2 3 0 4 554 5 0 x2 x1
3
32
2
6
22
22
4
54
2 3 0 4 5 5 4 5: 0 x_2 x_1
3
32
x1 x2
3
2
554
0
3 5;
10 (A.45)
586
Appendix: Numerical methods solving finite element dynamic equations
FIGURE A.1 The flowchart of the PTEM. PTEM, Program of Time Element Methods.
Appendix: Numerical methods solving finite element dynamic equations
587
TABLE A.1 Displacements calculated in Section A.7.1. t Δt
Newmark
Wilson-θ
Hilber-α
Hilber collocation
Central difference
Theoretical
0
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
2
1.4141
1.4256
1.4139
1.4142
1.4283
1.4142
4
1.0195
1.0532
1.0196
1.0196
1.0095
1.0008
6
0.0419
0.0974
0.0426
0.0420
20.0080
0.0017
8
20.9598
20.9063
20.9586
20.9596
201.0207
20.9984
10
21.4125
21.3996
21.4119
21.4125
21.4281
21.4142
12
21.0574
21.1115
21.0582
21.0574
20.9886
21.0024
TABLE A.2 Displacement x1 calculated in Section A.7.2. t Δt
Newmark
Wilson-θ
Hilber-α
Hilber collocation
Central difference
Theoretical
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
2
0.0504
0.0525
0.0507
0.0505
0.0307
0.0381
4
0.4846
0.4896
0.4847
0.4845
0.4871
0.4860
6
1.5805
1.5425
1.5790
1.5803
1.7009
1.6570
8
2.7607
2.6702
2.7578
2.7605
2.9134
2.8068
10
2.8505
2.8182
2.8502
2.8506
2.7711
2.8057
11
2.2840
2.3340
2.2865
2.2844
2.0368
2.1306
The theoretical solution of this problem is given by
1 5 x1 5 x2 3 5
22 4
"
pffiffiffi # 1 2 cos 2t pffiffiffi ; 1 2 cos 5t
(A.46)
pffiffiffi pffiffiffi which is a summation of two vibrations of periods T1 5 2π and T2 5 2π= 5. The time step Δt 5 0.1T2 is chosen in the calculations. The results obtained using PTEM are listed in Table A.2. The solutions for the Newmark and Wilson-θ methods are exactly same as the ones reported by Bathe and Wilson (1976) and by Bathe (1982, 1986).
A.7.3
Example 3
Here, a weak nonlinear problem of a system with 100 DOFs, as shown in Fig. A.2, is considered. For the ith degree of the system, the mass mi and the spring force fi are given as mi 5 1; fi 5 δi 1 0:2δ3i ;
(A.47)
where δi represents the elongation of the ith spring. Therefore the spring force of each spring is nonlinear. The Rayleigh damping of the system is assumed, given by C 5 0:01M 1 0:005K:
(A.48)
588
Appendix: Numerical methods solving finite element dynamic equations
C1,K1,m1
C100,K100,m100
•••
•
FIGURE A.2 A nonlinear system of 100 degrees of freedom.
TABLE A.3 Displacement x100 calculated in Section A.7.3. t Δt 0 48
Newmark
Wilson-θ
Hilber-α
20.000
20.000
20.000
8.7725
8.7537
Hilber collocation 20.000
8.7732
8.7727
Central difference 20.000 8.7838
Zheng and Tan 20.000 8.7210
108
23.130
23.245
23.132
23.131
23.028
23.123
149
4.008
3.879
4.006
4.007
4.113
4.080
206
23.597
23.588
23.596
23.596
23.646
23.622
247
2 6.765
26.740
26.765
26.766
26.857
26.862
287
26.985
26.928
26.984
26.985
27.088
26.988
The dynamic response of the system excited by five forces πt ; i 5 20; 40; 60; 80; 100 Ri 5 0:005i sin 50
(A.49)
is calculated using PTEM. The initial conditions are xi ð0Þ 5 0:2i;
x_i ð0Þ 5 0;
ði 5 1; 2; 3; . . .; 100Þ :
(A.50)
If the nonlinear part of total stiffness is relatively smaller than its linear part, we consider a constant linear stiffness 2 3 2 21 6 21 & & 7 7 K56 (A.51) 4 & 2 215 21 1 s and a nonlinear force term F~ , in which an element is given by h i s F~ i 5 0:2 ðxi 2xi21 Þ3 1 ðxi 2xi11 Þ3 ;
x0 5 x101 5 0:
(A.52)
The equation describing the system is now written as s M€x 1 C_x 1 Kx 5 R 1 F~ :
(A.53)
The time steps Δt 5 0.25 for the central difference method and Δt 5 2.5 for the other four methods are used. The parameters for each time integration method are the same as in Sections A.7.1 and A.7.2. The results obtained are shown and compared in Table A.3, in which the results by Zheng and Tan (1985) and Zheng (1986) are listed.
Appendix: Numerical methods solving finite element dynamic equations
A.8
Program of Time Element Methods
589
590
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
591
592
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
593
594
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
595
596
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
597
598
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
599
600
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
601
602
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
603
604
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
605
606
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
607
608
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
609
610
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
611
612
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
613
614
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
615
616
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
A.9
617
User manual
The variables used in the input file INPUT.DAT follow the rules of FORTAN. The users need to generate this INPUT. DAT as their input file to use the program. In Section A.10, the input files for the three examples are listed. The structure of the input file is defined as follows: Line 1: Maximum 72 letters allowed for users to define the title of problem Line 2: Main control data of 20 integer variables IJK, JC, NN, NWM, NWC, NWK, NITER, NPRT, NINT, NC, NO, NLD, NLDT, NLDU, NLDV, NAD, NAT, JCM, JCC, JCK The physical meanings of these variables are: IJK: Time element method control 5 1 Newmark method 5 2 Wilson-θ method 5 3 Hilber-α modification method 5 4 Hilber allocation method 5 5 Central difference method
618
Appendix: Numerical methods solving finite element dynamic equations
JC: Solution method control with the definition 51 M, C, K matrices are not changed at each time step, and loadR(t) is time function only, no iteration for the solution 52 M, C, K matrices are not changed at each time step, but loadR depends on the unknown solution and iteration needed except using the central difference method 53 M, C, K matrices are changed, which is required to be input at each time step, the load R(t) is time function only and no iteration needed NN: NWM: NWC: NWK: NITER: NPRT:
The order of matrices M, C, and K The length of matrix M in a one-dimensional array form The length of matrix C in a one-dimensional array form The length of matrix K in a one-dimensional array form Maximum iteration time number, automatically set to 30, if it is 0 Print control 5 1 Print the results during the process of calculation 5 2 Not print the results during the process of calculation NNIT: Load function control 5 0 Load not involving acceleration 5 2 Load involving acceleration (option not for this version) NC: Time step number for output results NO: Number of DOF that the results require NLD: Number of load functions, NLD 5 NLDT 1 NLDU 1 NLDV NLDT: Number of load functions as time function only NLDU: Number of nonlinear spring forces in the form fI 5 AI ðxI 2xE Þ2 1 BI ðxI 2xE Þ3 ; AI ; BI : constants to be given by users; xI : displacement at end I of spring; where force acted; xE : displacement at another end of sprint:
NLDV: Number of nonlinear damping forces in the form fI 5 AI ðx_I 2 x_E Þ2 1 BI ðx_I 2 x_E Þ3 ; AI ; BI : constants to be given by users; x_I : velocity at end I of damper; where damping force acted; x_E : velocity at another end of damper: NAD: Number of DOF where loads acted. If a DOF has several different forces, the total number of different forces minus 1 should be added into NAD NAT: Number of loads acting times JCM: M matrix change control at each step 5 1 Change 5 0 Not change JCC: C matrix change control at each step 5 1 Change 5 0 Not change JCK: K matrix change control at each step 5 1 Change 5 0 Not change
Appendix: Numerical methods solving finite element dynamic equations
619
Line 3: Eight parameters of time elements TS, TE, DT, GAMMA, BETA, THETA, ALPHA, TOL, defined as: TS: The starting time of solution TE: The ending time of solution DT: Time step (second) GAMMA: Parameter γ BETA: Parameter β THETA: Parameter θ ALPHA: Parameter α TOL: Convergence accuracy, which is set to 1025, if given as 0 Line 4: Loads arriving time data, total NAT real numbers for array AT(NAT) Time load function cards: Total NLDT sets of data. If NLDT 5 0, it is not needed. For each time function, the following data lines are required. Line 1:
Line 2: Line 3:
Line 4:
NLP: Load function control 5 0 Sinusoidal function .0 Discrete time function, NLP equals the data number SFTR: A coefficient factor Title of function: Total maximum 60 letters Sinusoidal function (only for NLP 5 0 case) in the form f ðtÞ 5 SFTR 3 SINð2π 3 FREQ 3 t 1 2π 3 PHASEÞ FREQ: Frequency (Hz) PHASE: Phase angle 5 phase degree/360 Discrete time function data (only for NLP . 0); total NLP sets data of time T and function value for the array (T(J), P(J), J 5 1, NLP), of which the ending time T(NLP) . TE, the ending time of solution
Nonlinear spring force data cards (only used when NLDU . 0): Total NLDU sets data into array AU(NLDU, 2) in the form ðAUðI; 1Þ; AUðI; 2Þ; I 5 1; NLDUÞ; AUðI; 1Þ 5 AI ; AUðI; 2Þ 5 BI : Nonlinear damping force data cards (only used when NLDV . 0): Total NLDV sets data into array AV(NLDV, 2) in the form ðAVðI; 1Þ; AVðI; 2Þ; I 5 1; NLDVÞ; AVðI; 1Þ 5 AI ; AVðI; 2Þ 5 BI : Original matrices data cards: 1.
Mass matrix RM(NWM) and its address array IM(NNN), NNN 5 NN 1 1, NWM 5 IM (NNN) 2 1 RM(NWM): One-dimensional stored mass matrix IM(NNN): The address of diagonal elements of mass matrix M in RM array
2.
Damping matrix RC(NWC) and its address array IC(NNN) (only used for NWC . 0 case) NNN 5 NN 1 1, NWC 5 IC(NNN) 2 1 RC(NWC): One-dimensional stored damping matrix IC(NNN): The address of diagonal elements of damping matrix C in RC array
620
Appendix: Numerical methods solving finite element dynamic equations
3.
Stiffness matrix RK(NWK) and its address array IK(NNN), NNN 5 NN 1 1, NWK 5 IK(NNN) 2 1 RK(NWK): One-dimensional stored stiffness matrix IK(NNN): The address of diagonal elements of stiffness matrix K in RK array
To illustrate this data structure, we give a symmetric matrix M and its corresponding arrays RM and IM as follows: 2 3 a e 0 0 6 b 0 07 7; RM 5 ða; b; c; d; e; f Þ; IM 5 ð1; 2; 4; 5; 7Þ; NNN 5 5: M56 4 d r5 e Here, for example, the third diagonal element d of M matrix locates at the fourth position in array RM, so that IM (3) 5 4. IM(5) 5 Length of RM 1 1 5 7. Initial conditions cards: U0(NN): Initial displacement vector, NN number of real numbers U1(NN): Initial velocity vector, NN umber of real numbers Load information data matrix NADF (NAD,4) and factor PFCTOR (NAD): Total NAD lines of which each line inputs the following five data elements: NADF(I,1): NADF(I,2):
NADF(I,3): NADF(I,4): PFCTOR(I):
Degree number at which Load I acts. This can be repeated to deal with the case if the same degree number is subject to different loads. Load I’s function ID number, which is in the order: NLDT time functions,NLDU nonlinear spring forces, NLDV nonlinear damping forces. Starting time number of load I, which is an integer of 1BNAT For the time function, set to 0; for nonlinear force function I at degree NADF(I,1), here input another end degree number. Load I’s coefficient factor
Elements change data of matrices RM, RC, RK: From the second time step to the NT step, calculated by NT 5 IDINT ((TE 2 TS)/DT) 2 1 NOUT 5 NT/NC 1 1 NT 5 (NOUT 2 1) 3 NC IF (IJK.EQ.5) NT 5 NT 1 1 to input the following data: JCM 5 1 input RM (NWM) JCC 5 1 input RC (NWC) JCK 5 1 input RK (NWK)
A.10
Input and output files for Examples 1 and 2 in Section A.7
As examples, we give the input and output files for Examples 1 and 2 discussed in Section A.7, which should help readers use the program manual for their problems of interest. The data files for Example 3 are quite a bit longer, so they are left out here.
Appendix: Numerical methods solving finite element dynamic equations
A.10.1
Example 1
INPUT.DAT Input file for using the Newmark method
621
622
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
623
624
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
625
626
Appendix: Numerical methods solving finite element dynamic equations
Numerical methods solving finite element dynamic equations This appendix gives some numerical methods, developed by Xing (2005), in FORTRAN designed to solve dynamic finite element (FE) equations. The principles, solution processes, examples, and comparison of five time integration methods—Newmark, Wilson-θ, Hilber-α, Hilber collocation, and central difference—which are often used to solve dynamic FE equations, are presented. A computer code Program of Time Element Methods (PTEM), developed by the author, is introduced and presented. This code provides a generalized functional module using these five algorithms to solve the FE equation t M x€ 1 t C x_ 1 t Kx 5 Rðt; x; x_ Þ describing linear and nonlinear structural dynamics. This functional program module is conveniently incorporated into any computer program used to investigate complex dynamics problems in maritime engineering.
A.1
Introduction
Direct integration methods, or time element methods, are important and necessary approaches in analyzing dynamic response problems of linear and nonlinear dynamic systems (Bathe and Wilson, 1976; Bathe, 1982, 1986, 1996; Belytschko and Hughes 1983; Zheng and Tan, 1985; Zheng 1986). Essentially, direct numerical integration methods are based on two ideas. First, instead of satisfying the FE dynamic equation at any time t, these methods are aimed to satisfy it only at discrete time intervals Δt apart. The second is that a variation of displacements, velocities, and accelerations within each time interval is assumed. In depending on the different discrete points and variation assumptions used, there are different time integration schemes with different performances of accuracy, stability, and cost of solution. The central difference method is an often used explicit method with three-point interpolation. Houbolt (1950) developed a four-point implicit scheme. Newmark (1959) introduced two parameters γ and β to represent displacement and velocity, respectively, and created a very useful algorithm called the Newmark method. Hilber and Taylor (1977) introduced another parameter, α, in order to control numerical damping and to improve accuracy, into the Newmark scheme and created the Hilber-α method. Wilson et al. (1973), assuming a linear acceleration distribution over a time interval (t,t 1 θΔt), constructed the Wilson-θ method. Combining the ideas of the Wilson and Newmark methods, Hilber and Hughes (1978) developed another method called the Hilber collocation method. Zienkiewiez (1977) investigated the constructions of the Newmark and other time integration methods using a weighted residual approach. Based on variational principles and interpolations in the time domain, Xing and Zheng (1985, 1987) studied existing methods and created a new approach. In an integrated variational frame for the interpolations in space and time domains, the concept of a four-dimensional FE method was proposed. In the current FE software, the five methods—the central difference, Newmark, Wilson-θ, Hilber-α, Hilber collocation schemes—are widely used. This appendix gives a description and comparison of these five methods. A computer code PTEM, developed by the author, is introduced and listed here for application by readers.
A.2 A.2.1
Fundamental iteration formulations Generalized dynamic equation
In a linear FE analysis (FEA), a generalized dynamic equation is expressed as M€x 1 C_x 1 Kx 5 R;
(A.1)
577
578
Appendix: Numerical methods solving finite element dynamic equations
where M, C, and K are constant mass, damping, and stiffness matrices of the system, x, x_ , and x€ , represent the displacement, velocity, and acceleration vectors, and R is a time-dependent load vector. For nonlinear systems, Eq. (A.1) can be replaced by the following increment form: t
MΔ€x 1 t CΔ_x 1 t KΔx 5 Δt R;
(A.2)
where tM, tC, and tK denote the tangent mass, damping, and stiffness matrices at time t, and Δt R is the load increment. Since at each time step, the tangent mass, damping, and stiffness matrices are recalculated, the cost of solving Eq. (A.2) is quite high. For practical problems, nonlinearity is mainly caused by stiffness and damping but not involving the mass of the system. To reduce the computational cost, Eq. (A.2) is represented as t 1 Δt
M x€ ðkÞ 1 0 CΔ_xðkÞ 1 0 KΔxðkÞ 5 t 1 Δt R 2 t 1 Δt Fd ð_xðk21Þ Þ 2 t 1 Δt Fs ðxðk21Þ Þ;
t 1 Δt ðkÞ
5 t 1 Δt xðk 2 1Þ 1 ΔxðkÞ ;
t 1 Δt ðkÞ
5 t 1 Δt _x ðk 2 1Þ 1 Δ_xðkÞ ;
x
x_
(A.3)
k 5 1; 2; 3; :::
where 0C and 0K are the initial damping and stiffness matrices at time t 5 0, and Fd and Fs represent the damping and stiffness forces that can be expressed in a summation of linear and nonlinear components as follows: d Fd 5 0 C x_ 1 F~ ð_xÞ;
s Fs 5 0 Kx 1 F~ ðxÞ:
(A.4)
Eq. (A.3) is now rewritten in the form t 1 Δt
t 1 Δt ~ R 2 t 1 Δt F~ d ð_xðk21Þ Þ 2 t 1 Δt F~ s ðxðk21Þ Þ 5 t 1 Δt R; ðM€xðkÞ 1 C_xðkÞ 1 KxðkÞ Þ 5
(A.5)
where for convenience the superscripts of the damping and stiffness matrices are neglected. Eq. (A.5) provides an iteration formulation to obtain a solution x(k) at iteration k from previous result x(k21). The following three convergence criteria may be chosen to check whether the convergence is reached at iteration k.
A.2.2
Iteration criterions
A.2.2.1 Energy criterion ΔxðkÞ t 1 Δt ½R 2 Fd ð_xðk21Þ 2 Fs ðxðk21Þ 2 M€xðk21Þ Δxð1Þ U½t 1 Δt R 2 t Fd 2 t Fs 2 Mt1Δt x€ ð0Þ
(A.6)
# ETOL:
A.2.2.2 Force criterion :t1Δt ½R2Fd ð_xðk21Þ 2Fs ðxðk21Þ 2M€xðk21Þ :2 :½t1Δt R2 t Fd 2 t Fs 2Mt1Δt x€ ð0Þ :2
(A.7)
# RTOL:
A.2.2.3 Displacement criterion :ΔxðkÞ :2 # DTOL: :xðkÞ :2
Here :a:2 5
pffiffiffiffiffiffiffiffi aUa, ETOL, RTOL, and DTOL represent accuracy errors defined by users.
(A.8)
Appendix: Numerical methods solving finite element dynamic equations
A.3 A.3.1
579
Five numerical integration schemes Newmark method
The Newmark method assumes that the velocity and displacement vectors at time step n 1 1 take the forms x_ n11 5 x_ n 1 ½ð1 2 γÞ€xn 1 γ x€ n11 ; xn11 5 xn 1 x_ n Δt 1 ½ð0:5 2 βÞ€xn 1 β x€ n11 Δt2 ;
(A.9)
~ n11 ; M€xn11 1 C_xn11 1 Kxn11 5 R from which the displacement, velocity, and acceleration vectors xn11, x_ n11 , and x€ n11 at time step n 1 1 are derived from their values at the previous time step n. The unconditional stable condition of the Newmark method requires that β $ 0:25ð0:51γÞ2 :
γ $ 0:5;
A.3.2
(A.10)
Wilson-θ method
This method assumes that the acceleration in time interval τAðt; t 1 θΔtÞ varies in a linear form: τ ð€xn1θΔt 2 x€ n Þ: x€ n1τ 5 x€ n 1 θΔt
(A.11)
The integration of this equation gives that x_ n1τ 5 x_ n 1 x€ n τ 1
τ2 ð€xn1θΔt 2 x€ n Þ; 2θΔt
τ2 τ3 xn1τ 5 xn 1 x_ n τ 1 x€ n 1 ð€xn1θΔt 2 x€ n Þ; 2 6θΔt
(A.12)
which, at τ 5 θΔt, x_ n1θΔt , and x€ n1θΔt , can be solved in terms of xn1θΔt . The substitution of the results of x_ n1θΔt and x€ n1θΔt into the dynamic equation M€xn1θΔt 1 C_xn1θΔt 1 Kxn1θΔt ~ n1θΔt 5 R ~ n 1 θðR ~ n1Δt 2 R ~ nÞ 5R
(A.13)
provides a solution for xn1θΔt . The solution of xn1Δt , x_ n1Δt , and x€ n1Δt can be obtained by 3 6 6 x€ n 2 2 x_ n 1 3 2 ðxn1θΔt 2 xn Þ; x€ n11 5 1 2 θ θ Δt θ Δt x_ n11 5 x_ n 1
(A.14a)
Δt ð€xn11 1 x€ n Þ; 2
Δt2 ð€xn11 1 2€xn Þ: xn11 5 xn 1 Δtx_ n 1 6
(A.14b)
The unconditional stable condition of the method requires θ $ 1:37.
A.3.3
Hilber-α method
Based on the Newmark method, Hilber introduced a parameter α into Eq. (A.9) and obtained that ~ n11 ; M€xn11 1 C_xn11 1 K½xn11 ð1 1 αÞ 2 αxn 5 R
(A.15)
580
Appendix: Numerical methods solving finite element dynamic equations
developing a modification algorithm. The parameter α controls human-made damping and increases the accuracy. The unconditional stable condition of this method requires that γ 5 0:5 2 α; β 5 0:25ð1 2 αÞ;
A.3.4
20:5 # α # 0 :
(A.16)
Hilber collocation method
Combining the Wilson-θ and Newmark methods, Hilber constructed a collocation scheme using the following form: x€ n1θ 5 x€ n 1 θð€xn11 2 x€ n Þ; x_ n1θ 5 x_ n 1 θΔt½ð1 2 γÞ€xn 1 γ x€ n1θ ; xn1θ 5 xn 1 θΔtx_ n 1 ðθΔtÞ2 ½ð0:5 2 βÞ€xn 1 β x€ n1θ ; ~ n1θ ; M€xn1θ 5 C€xn1θ 1 Kxn1θ 5 R
(A.17)
in which x_ n1θ and x€ n1θ are represented in terms of xn1θ , and then xn1θ is obtained by solving ^ ^ n1θ 5 R; Kx 1 γ ^5 C 1 K; M1 K βθΔt βðθΔtÞ2 ^ 5 ð1 2 θÞR ~ n 1 θR ~ n11 ; R 0 1 3 2 1 1 1 xn 1 1 M4 x_ n 1 @ 2 1Ax€ n 5; βðθΔtÞ 2β βðθΔtÞ2 0 1 0 1 2 3 γ γ γ xn 1 @ 2 1Ax_ n 1 @ 2 1AθΔtx€ n 5: 1 C4 βðθΔtÞ β 2β
(A.18)
Hilber demonstrated that this method is unconditionally stable provided that γ 5 0;
θ 2θ2 2 1 $β$ ; 2ðθ 1 1Þ 4ð2θ3 2 1Þ
θ$1
(A.19)
This method reduces to a Newmark method of γ 5 0:5 and β 5 0:25 if θ 5 1 and to a Wilson-θ method if γ 5 1=2 and β 5 1=6.
A.3.5
Central difference method
For the central difference method, it is assumed that
x_ n 5
1 ðxn11 2 xn21 Þ; 2Δt (A.20a)
1 x€ n 5 2 ðxn11 2 2xn 1 xn21 Þ; Δt ~ n; M€xn 1 C_xn 1 Kxn 5 R
(A.20b)
Appendix: Numerical methods solving finite element dynamic equations
581
from which it follows that ^ n11 5 R; ^ Kx ^ 5 1 M 1 1 C; K Δt2 2Δt 0 1 0 1 2 1 1 ^ 5R ~ n 2 @K 2 CAxn21 ; R M A xn 2 @ 2 M 2 Δt2 Δt 2Δt x21 5 x0 2 Δtx_ 0 1
(A.21)
Δt2 x€ 0 : 2
The central difference method is an explicit method for which iteration is not needed, and therefore it is widely used in nonlinear analysis. This method is conditionally stable, and the stable condition requires Δt #
Tm ; π
(A.22)
where Tm is the minimum vibration period of the problem.
A.4
Stability and accuracy
As previously described, the four methods, except for the central difference method, are unconditionally stable, implying that for any time step Δt, an approximate solution can be definitely obtained. However, for a large time step Δt, the accuracy of the result is lower but at less of a cost of calculation. Here, the stability and accuracy analysis of the methods are given for a reference in applications.
A.4.1
Stability and spectral radius
The stability of an integration method is determined by examining the behavior of the numerical solution for arbitrary initial conditions. Assume that the iteration formulation of a method to determine the solution t 1 Δt x^ from the solution t x^ is defined as t 1 Δt
x^ 5 At x^ ;
(A.23)
where A represents the matrix operator of a method. Let ρ(A) be the spectral radius of A defined as ρðAÞ 5 max jλj j;
j 5 1; 2; 3; . . . :
(A.24)
The necessary and sufficient condition of a stable time integration method is ρðAÞ # 1:
A.4.2
(A.25)
Period elongation and amplitude decay
Various studies of the accuracies of time integration methods have been reported. Here, we only summarize some important solution characteristics. For this purpose, we consider the solution of the initial value problem defined by x€ 1 ω2 x 5 0;
0
x 5 1:0; 0 x_ 5 0:0;
(A.26)
for which the exact solution is x 5 cos ωt. However, there is a period elongation PE and an amplitude decay AD in the approximate solution obtained by time integration methods. It has been defined that Percentage amplitude decay 5 AD 3 100%; Percentage period elongation 5
PE 3 100%; T
(A.27)
582
Appendix: Numerical methods solving finite element dynamic equations
to describe the accuracy of a method. Also, a damping factor η has been introduced as AD 5 1 2 e22πη ;
η5
2lnð1 2 ADÞ : 2π
(A.28)
For the often used time integration methods, the references by Hilber et al. (1977) and Hilber and Hughes (1978) present the spectral radii ρ(A), the percentage amplitude decay, the damping factor, and the percentage of period elongation as functions of Δt/T. These curves provide a good comparison of the discussed time integration methods. From the point of view of engineering applications, all of the five methods can provide sufficient accuracy for engineering analysis, provided a suitable time step chosen.
A.5
Implementation
The implementation process of these five time integration methods is described as follows.
A.5.1
Initial calculations
1. Form mass, stiffness, and damping matrices M, K, and C. 2. Initialize 0 x, 0 x_ , and 0 x€ 5 M21 ð0 R 2 C0 x_ 2 K0 xÞ. For the central difference method, you need to calculate 21 x 5 0 x 5 0 x_ Δt 1 0:50 x€ Δt2 . 3. Select time step Δt , Δtcr , and calculate the integration constants.
A.5.1.1 Newmark method γ $ 0:5; Ak 5 1;
β $ 0:25ð0:51γÞ2 ; A0 5 1ðβΔt2 Þ; A1 5 γ=ðβΔtÞ;
A2 5 1=ðβΔtÞ;
A3 5 1ð2βÞ 2 1;
A5 5 ðγ=β 2 2ÞΔt=2;
A4 5 γ=β 2 1;
A6 5 ð1 2 γÞΔt;
A7 5 γΔt:
A.5.1.2 Hilber-α method 0:5 # α # 0;
γ 5 0:5 2 α; β 5 0:25ð12αÞ2 ;
Ak 5 1 1 α;
A0 5 1=ðβΔt2 Þ;
A2 5 1=ðβΔtÞ;
A1 5 γ=ðβΔtÞ;
A3 5 1=ð2βÞ 2 1;
A5 5 ðγ=β 2 2ÞΔt=2;
A4 5 γ=β 2 1;
A6 5 ð1 2 γÞΔt;
A7 5 γΔt:
A.5.1.3 Wilson-θ method θ $ 1:37;
Ak 5 1;
A0 5 6=ðθΔt2 Þ; A3 5 θΔt=2; A6 5 1 2 3=θ;
A1 5 3=ðθΔtÞ; A2 5 2A1 ;
A4 5 A0 =θ;
A5 5 2 A2 =θ;
A7 5 Δt=2; A8 5 Δt2 =6:
Appendix: Numerical methods solving finite element dynamic equations
583
A.5.1.4 Hilber collocation method θ $ 1;
γ 5 0:5;
Ak 5 1;
θ 2θ2 2 1 $β$ ; 2ðθ 1 1Þ 4ð2θ3 2 1Þ
A0 5 1=ðβθ2 Δt2 Þ;
A1 5 γ=ðβθΔtÞ;
A2 5 1=ðβθΔtÞ; A3 5 1=ð2βÞ 2 1; A4 5 γ=β 2 1; γ A0 21 2 1 θΔt; A6 5 A5 5 ; A7 5 2 ; 2β θ βθ Δt A8 5 2
1 ; 2βθ
B1 5 Δtγ;
A9 5 ð1 2 γÞΔt;
B2 5 ð0:5 2 βÞθ2 Δt2 ;
B3 5 βθ2 Δt2 :
A.5.1.5 Central difference method Δt # Tn =π 5 Δtcr ; Ak 5 0; A0 5 1=Δt2 ;
A.5.2
A1 5 1ð2ΔtÞ;
A2 5 2A0 :
Form effective stiffness matrix ^ 5 A0 M 1 A1 C 1 Ak K: K
A.5.3
Triangular decomposition of matrix K^ ^ 5 LDLT : K
A.5.4
(A.29a)
(A.29b)
Calculations at each time step
1. Form effective loads at time t 1 Δt. For a convenience in equation expression, we define that t
R 5 MðA0 t x 1 A2 t x_ 1 A3 t x€ Þ 1 CðA1 t x 1 A4 t x_ 1 A5 t x€ Þ:
(A.30)
Newmark method t 1 Δt
^ 5 t 1 Δt R 1 t R: R
(A.31)
t 1 Δt
^ 5 t 1 Δt R 1 t R: R
(A.32)
Hilber-α method
Wilson-θ method t 1 Δt
^ 5 A2 t R 1 θt 1 Δt R 2 A2 Mt x€ R 2 Kðt x 1 A3 t x_ 1 A5 t x€ Þ 2 Cðt x 1 A4 t x€ Þ:
(A.33)
Hilber collocation method t 1 Δt
^ 5 ð1 2 θÞt R 1 θt 1 Δt R 1 t R: R
(A.34)
584
Appendix: Numerical methods solving finite element dynamic equations
Central difference method t 1 Δt
^ 5 t R 1 MðA2 t x 2 A0 t 2 Δt xÞ 2 Kt x 1 CA1 t 1 Δt x: R
(A.35)
2. Solve for x^ at time t 1 Δt: ^ ^ x 5 LDLT x^ 5 t 1 Δt R: K^
(A.36)
3. Evaluate displacements, velocities, and accelerations at time t 1 Δt. Newmark and Hilber-α methods t 1 Δt
x 5 x^ ;
t 1 Δt
x€ 5 A0 ð^x 2 t x^ Þ 2 A2 t x_ 2 A3 t x€ ;
(A.37)
x 5 t x_ 1 A6 t x€ 1 A7 t 1 Δt x€ :
t 1 Δt _
Wilson-θ method t 1 Δt
x€ 5 A4 ð^x 2 t x^ Þ 1 A5 t x_ 1 A6 t x€ ; x 5 t x_ 1 A7 ðt x€ 1 t 1 Δt x€ Þ;
t 1 Δt _ t 1 Δt
(A.38)
x 5 t x 1 Δtt x_ 1 A8 ð2t x€ 1 t 1 Δt x€ Þ:
Hilber collocation method t 1 Δt
x€ 5 A6 ð^x 2 t x^ Þ 1 A7 t x_ 1 A8 t x€ ; x 5 t x_ 1 A9 t x€ 1 B1 t 1 Δt x€ ;
t 1 Δt _ t 1 Δt
x 5 x 1 Δt x 1 B2 x€ 1 B3 t
t_
t
(A.39) t 1 Δt
x€ :
Central difference method t 1 Δt
x 5 x^ ;
x 5 A1 ðt 1 Δt x 2 t 2 Δt xÞ;
t_ t
A.6
(A.40)
x€ 5 A0 ðt 1 Δt x 2 2t x 1 t 2 Δt xÞ:
Time element program
The time element program PTEM, developed by Xing (1989), is a functional module that includes the five time integration methods just discussed. This program is used to solve linear and nonlinear dynamic FEA equations to obtain displacements, velocities, and accelerations of the system at time t 1 Δt from their values at time t. It can be incorporated into an FEA program as a functional module or directly used to solve dynamic equations obtained using other software.
A.6.1
Main characteristics
1. The program consists of about 1200 FORTRAN sentences. All functions of the program are completed in the module. It provides a convenient interface for users to connect the module into other FEA programs. 2. There are three forms to solve Eq. (A.1), which is controlled by a control variable JC as follows: JC 5 1: The matrices M, C, and K are assumed unchanged at each time step, and load R is dependent on time only. This function is mainly used to solve linear problems with no iteration required at each time step. JC 5 2: The matrices M, C, and K are still assumed unchanged at each time step, but the load R is dependent on unknown displacement and velocity. For the four time integration methods except the central difference method, iteration is required at each time step. Users can choose the convergence criteria. This process is mainly used to solve problems with weak nonlinearity. For problems with strong nonlinearity, iteration convergence may not reach, and the central difference method can be used.
Appendix: Numerical methods solving finite element dynamic equations
585
JC 5 3: The matrices M, C, and K are changed at each time step, but the load R is independent of only the time function unknown variables. During the calculation at each time step the matrices M, C, and K are required to be updated using FEA software. There are three flags JCM, JCC, and JCK to control the changes of matrices M, C, and K, respectively. 3. A control flag IJK is designed for users to choose one of the five time integration methods. 4. The program can automatically check the parameters α, β, γ, and θ to satisfy the stability conditions. 5. The load functions designed in the program include sinusoidal functions, discrete time functions given by input data, and the nonlinear spring function FUNC1 5 Ax2 1 Bx3
(A.41)
FUNC1 5 A_x2 1 B_x3 ;
(A.42)
and nonlinear damping function
where A and B are constants defined by users. Users can input their load functions, defined by themselves, and define the starting time of each load. 6. Output results are the displacements, velocities, and accelerations at the selected degrees of freedom (DOFs) by users. Users can choose the time intervals of the output data.
A.6.2
Flowchart of the program
A flowchart of the program is shown in Fig. A.1. All the input data, the values of control variables, and the method parameters are supplied by an INPUT.DAT file prepared by users. A file OUTPUT.DAT produced automatically by the program gives the output results.
A.7 A.7.1
Examples Example 1
A nondamping free vibration of a system with 1 DOF is defined by the following equation: x€ 1 x 5 0;
xð0Þ 5 1:0;
_ 5 1:0 ; xð0Þ
(A.43)
which has a theoretical solution x5
pffiffiffi π 2 cos t 5 : 4
(A.44)
Table A.1 lists the displacements calculated using the program in association with the parameters Δt 5 π=8, γ 5 0:5, β 5 0:25, α 5 20:01, and θ 5 1:4 for the Wilson-θ method and θ 5 1:001 for the Hilber collocation method.
A.7.2
Example 2
Let us consider a system with 2 DOFs described by the following equation: 2 4
2 0 0 1
2
54
x€1 x€2
3
2
514
2 3 0 4 554 5 0 x2 x1
3
32
2
6
22
22
4
54
2 3 0 4 5 5 4 5: 0 x_2 x_1
3
32
x1 x2
3
2
554
0
3 5;
10 (A.45)
586
Appendix: Numerical methods solving finite element dynamic equations
FIGURE A.1 The flowchart of the PTEM. PTEM, Program of Time Element Methods.
Appendix: Numerical methods solving finite element dynamic equations
587
TABLE A.1 Displacements calculated in Section A.7.1. t Δt
Newmark
Wilson-θ
Hilber-α
Hilber collocation
Central difference
Theoretical
0
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
2
1.4141
1.4256
1.4139
1.4142
1.4283
1.4142
4
1.0195
1.0532
1.0196
1.0196
1.0095
1.0008
6
0.0419
0.0974
0.0426
0.0420
20.0080
0.0017
8
20.9598
20.9063
20.9586
20.9596
201.0207
20.9984
10
21.4125
21.3996
21.4119
21.4125
21.4281
21.4142
12
21.0574
21.1115
21.0582
21.0574
20.9886
21.0024
TABLE A.2 Displacement x1 calculated in Section A.7.2. t Δt
Newmark
Wilson-θ
Hilber-α
Hilber collocation
Central difference
Theoretical
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
2
0.0504
0.0525
0.0507
0.0505
0.0307
0.0381
4
0.4846
0.4896
0.4847
0.4845
0.4871
0.4860
6
1.5805
1.5425
1.5790
1.5803
1.7009
1.6570
8
2.7607
2.6702
2.7578
2.7605
2.9134
2.8068
10
2.8505
2.8182
2.8502
2.8506
2.7711
2.8057
11
2.2840
2.3340
2.2865
2.2844
2.0368
2.1306
The theoretical solution of this problem is given by
1 5 x1 5 x2 3 5
22 4
"
pffiffiffi # 1 2 cos 2t pffiffiffi ; 1 2 cos 5t
(A.46)
pffiffiffi pffiffiffi which is a summation of two vibrations of periods T1 5 2π and T2 5 2π= 5. The time step Δt 5 0.1T2 is chosen in the calculations. The results obtained using PTEM are listed in Table A.2. The solutions for the Newmark and Wilson-θ methods are exactly same as the ones reported by Bathe and Wilson (1976) and by Bathe (1982, 1986).
A.7.3
Example 3
Here, a weak nonlinear problem of a system with 100 DOFs, as shown in Fig. A.2, is considered. For the ith degree of the system, the mass mi and the spring force fi are given as mi 5 1; fi 5 δi 1 0:2δ3i ;
(A.47)
where δi represents the elongation of the ith spring. Therefore the spring force of each spring is nonlinear. The Rayleigh damping of the system is assumed, given by C 5 0:01M 1 0:005K:
(A.48)
588
Appendix: Numerical methods solving finite element dynamic equations
C1,K1,m1
C100,K100,m100
•••
•
FIGURE A.2 A nonlinear system of 100 degrees of freedom.
TABLE A.3 Displacement x100 calculated in Section A.7.3. t Δt 0 48
Newmark
Wilson-θ
Hilber-α
20.000
20.000
20.000
8.7725
8.7537
Hilber collocation 20.000
8.7732
8.7727
Central difference 20.000 8.7838
Zheng and Tan 20.000 8.7210
108
23.130
23.245
23.132
23.131
23.028
23.123
149
4.008
3.879
4.006
4.007
4.113
4.080
206
23.597
23.588
23.596
23.596
23.646
23.622
247
2 6.765
26.740
26.765
26.766
26.857
26.862
287
26.985
26.928
26.984
26.985
27.088
26.988
The dynamic response of the system excited by five forces πt ; i 5 20; 40; 60; 80; 100 Ri 5 0:005i sin 50
(A.49)
is calculated using PTEM. The initial conditions are xi ð0Þ 5 0:2i;
x_i ð0Þ 5 0;
ði 5 1; 2; 3; . . .; 100Þ :
(A.50)
If the nonlinear part of total stiffness is relatively smaller than its linear part, we consider a constant linear stiffness 2 3 2 21 6 21 & & 7 7 K56 (A.51) 4 & 2 215 21 1 s and a nonlinear force term F~ , in which an element is given by h i s F~ i 5 0:2 ðxi 2xi21 Þ3 1 ðxi 2xi11 Þ3 ;
x0 5 x101 5 0:
(A.52)
The equation describing the system is now written as s M€x 1 C_x 1 Kx 5 R 1 F~ :
(A.53)
The time steps Δt 5 0.25 for the central difference method and Δt 5 2.5 for the other four methods are used. The parameters for each time integration method are the same as in Sections A.7.1 and A.7.2. The results obtained are shown and compared in Table A.3, in which the results by Zheng and Tan (1985) and Zheng (1986) are listed.
Appendix: Numerical methods solving finite element dynamic equations
A.8
Program of Time Element Methods
589
590
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
591
592
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
593
594
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
595
596
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
597
598
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
599
600
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
601
602
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
603
604
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
605
606
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
607
608
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
609
610
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
611
612
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
613
614
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
615
616
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
A.9
617
User manual
The variables used in the input file INPUT.DAT follow the rules of FORTAN. The users need to generate this INPUT. DAT as their input file to use the program. In Section A.10, the input files for the three examples are listed. The structure of the input file is defined as follows: Line 1: Maximum 72 letters allowed for users to define the title of problem Line 2: Main control data of 20 integer variables IJK, JC, NN, NWM, NWC, NWK, NITER, NPRT, NINT, NC, NO, NLD, NLDT, NLDU, NLDV, NAD, NAT, JCM, JCC, JCK The physical meanings of these variables are: IJK: Time element method control 5 1 Newmark method 5 2 Wilson-θ method 5 3 Hilber-α modification method 5 4 Hilber allocation method 5 5 Central difference method
618
Appendix: Numerical methods solving finite element dynamic equations
JC: Solution method control with the definition 51 M, C, K matrices are not changed at each time step, and loadR(t) is time function only, no iteration for the solution 52 M, C, K matrices are not changed at each time step, but loadR depends on the unknown solution and iteration needed except using the central difference method 53 M, C, K matrices are changed, which is required to be input at each time step, the load R(t) is time function only and no iteration needed NN: NWM: NWC: NWK: NITER: NPRT:
The order of matrices M, C, and K The length of matrix M in a one-dimensional array form The length of matrix C in a one-dimensional array form The length of matrix K in a one-dimensional array form Maximum iteration time number, automatically set to 30, if it is 0 Print control 5 1 Print the results during the process of calculation 5 2 Not print the results during the process of calculation NNIT: Load function control 5 0 Load not involving acceleration 5 2 Load involving acceleration (option not for this version) NC: Time step number for output results NO: Number of DOF that the results require NLD: Number of load functions, NLD 5 NLDT 1 NLDU 1 NLDV NLDT: Number of load functions as time function only NLDU: Number of nonlinear spring forces in the form fI 5 AI ðxI 2xE Þ2 1 BI ðxI 2xE Þ3 ; AI ; BI : constants to be given by users; xI : displacement at end I of spring; where force acted; xE : displacement at another end of sprint:
NLDV: Number of nonlinear damping forces in the form fI 5 AI ðx_I 2 x_E Þ2 1 BI ðx_I 2 x_E Þ3 ; AI ; BI : constants to be given by users; x_I : velocity at end I of damper; where damping force acted; x_E : velocity at another end of damper: NAD: Number of DOF where loads acted. If a DOF has several different forces, the total number of different forces minus 1 should be added into NAD NAT: Number of loads acting times JCM: M matrix change control at each step 5 1 Change 5 0 Not change JCC: C matrix change control at each step 5 1 Change 5 0 Not change JCK: K matrix change control at each step 5 1 Change 5 0 Not change
Appendix: Numerical methods solving finite element dynamic equations
619
Line 3: Eight parameters of time elements TS, TE, DT, GAMMA, BETA, THETA, ALPHA, TOL, defined as: TS: The starting time of solution TE: The ending time of solution DT: Time step (second) GAMMA: Parameter γ BETA: Parameter β THETA: Parameter θ ALPHA: Parameter α TOL: Convergence accuracy, which is set to 1025, if given as 0 Line 4: Loads arriving time data, total NAT real numbers for array AT(NAT) Time load function cards: Total NLDT sets of data. If NLDT 5 0, it is not needed. For each time function, the following data lines are required. Line 1:
Line 2: Line 3:
Line 4:
NLP: Load function control 5 0 Sinusoidal function .0 Discrete time function, NLP equals the data number SFTR: A coefficient factor Title of function: Total maximum 60 letters Sinusoidal function (only for NLP 5 0 case) in the form f ðtÞ 5 SFTR 3 SINð2π 3 FREQ 3 t 1 2π 3 PHASEÞ FREQ: Frequency (Hz) PHASE: Phase angle 5 phase degree/360 Discrete time function data (only for NLP . 0); total NLP sets data of time T and function value for the array (T(J), P(J), J 5 1, NLP), of which the ending time T(NLP) . TE, the ending time of solution
Nonlinear spring force data cards (only used when NLDU . 0): Total NLDU sets data into array AU(NLDU, 2) in the form ðAUðI; 1Þ; AUðI; 2Þ; I 5 1; NLDUÞ; AUðI; 1Þ 5 AI ; AUðI; 2Þ 5 BI : Nonlinear damping force data cards (only used when NLDV . 0): Total NLDV sets data into array AV(NLDV, 2) in the form ðAVðI; 1Þ; AVðI; 2Þ; I 5 1; NLDVÞ; AVðI; 1Þ 5 AI ; AVðI; 2Þ 5 BI : Original matrices data cards: 1.
Mass matrix RM(NWM) and its address array IM(NNN), NNN 5 NN 1 1, NWM 5 IM (NNN) 2 1 RM(NWM): One-dimensional stored mass matrix IM(NNN): The address of diagonal elements of mass matrix M in RM array
2.
Damping matrix RC(NWC) and its address array IC(NNN) (only used for NWC . 0 case) NNN 5 NN 1 1, NWC 5 IC(NNN) 2 1 RC(NWC): One-dimensional stored damping matrix IC(NNN): The address of diagonal elements of damping matrix C in RC array
620
Appendix: Numerical methods solving finite element dynamic equations
3.
Stiffness matrix RK(NWK) and its address array IK(NNN), NNN 5 NN 1 1, NWK 5 IK(NNN) 2 1 RK(NWK): One-dimensional stored stiffness matrix IK(NNN): The address of diagonal elements of stiffness matrix K in RK array
To illustrate this data structure, we give a symmetric matrix M and its corresponding arrays RM and IM as follows: 2 3 a e 0 0 6 b 0 07 7; RM 5 ða; b; c; d; e; f Þ; IM 5 ð1; 2; 4; 5; 7Þ; NNN 5 5: M56 4 d r5 e Here, for example, the third diagonal element d of M matrix locates at the fourth position in array RM, so that IM (3) 5 4. IM(5) 5 Length of RM 1 1 5 7. Initial conditions cards: U0(NN): Initial displacement vector, NN number of real numbers U1(NN): Initial velocity vector, NN umber of real numbers Load information data matrix NADF (NAD,4) and factor PFCTOR (NAD): Total NAD lines of which each line inputs the following five data elements: NADF(I,1): NADF(I,2):
NADF(I,3): NADF(I,4): PFCTOR(I):
Degree number at which Load I acts. This can be repeated to deal with the case if the same degree number is subject to different loads. Load I’s function ID number, which is in the order: NLDT time functions,NLDU nonlinear spring forces, NLDV nonlinear damping forces. Starting time number of load I, which is an integer of 1BNAT For the time function, set to 0; for nonlinear force function I at degree NADF(I,1), here input another end degree number. Load I’s coefficient factor
Elements change data of matrices RM, RC, RK: From the second time step to the NT step, calculated by NT 5 IDINT ((TE 2 TS)/DT) 2 1 NOUT 5 NT/NC 1 1 NT 5 (NOUT 2 1) 3 NC IF (IJK.EQ.5) NT 5 NT 1 1 to input the following data: JCM 5 1 input RM (NWM) JCC 5 1 input RC (NWC) JCK 5 1 input RK (NWK)
A.10
Input and output files for Examples 1 and 2 in Section A.7
As examples, we give the input and output files for Examples 1 and 2 discussed in Section A.7, which should help readers use the program manual for their problems of interest. The data files for Example 3 are quite a bit longer, so they are left out here.
Appendix: Numerical methods solving finite element dynamic equations
A.10.1
Example 1
INPUT.DAT Input file for using the Newmark method
621
622
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
623
624
Appendix: Numerical methods solving finite element dynamic equations
Appendix: Numerical methods solving finite element dynamic equations
625
626
Appendix: Numerical methods solving finite element dynamic equations