The following transient analysis topics are available:
The transient analysis solution method (ANTYPE,TRANS) used depends on the DOFs involved. Structural, acoustic, and other second order systems (that is, the systems are second order in time) are solved using one method and the thermal, magnetic, electrical and other first order systems are solved using another. Each method is described subsequently. If the analysis contains both first and second order DOFs (for example, structural and magnetic), then each DOF is solved using the appropriate method. For matrix coupling between first and second order effects such as for piezoelectric analysis, a combined procedure is used.
Initial conditions are known.
Gyroscopic or Coriolis effects are included in a structural analysis when requested (using the CORIOLIS command).
For most structural dynamics problems of a mechanical system, the spatial discretization for the principle of virtual work using the finite element method gives the finite element semi-discrete equation of motion as follows:
(15–5) |
where:
= structural mass matrix |
= structural damping matrix |
= nodal acceleration vector |
= nodal velocity vector |
= nodal displacement vector |
= internal load vector |
= applied load vector |
Three methods are available for solving Equation 15–5:
Newmark time integration method -- Used for implicit transient analyses as described below. This method is requested by setting
TINTOPT
= NMK (which is the default) on the TRNOPT command.HHT time integration method -- Used for implicit transient analyses as described below. This method is an extension of the Newmark time integration method and is requested by setting
TINTOPT
= HHT on the TRNOPT command.Backward Euler time integration method -- Used for implicit transient analyses as described below. This method can be requested by issuing the command TINTP,QUAS.
The structural dynamics problems concerned with the mechanical behavior governed by the above differential equation can be classified into two classes; that is, linear and nonlinear problems.
In linear structural dynamics systems, the internal load is linearly proportional to the nodal displacement, and the structural stiffness matrix remains constant. Therefore, Equation 15–5 can be rewritten as:
(15–6) |
where:
= structural stiffness matrix |
Among direct time integration methods for numerically solving the finite element semi-discrete equation of motion given in Equation 15–6, several methods such as the Newmark method (Newmark([402])) and the generalized- method (Chung and Hulbert([348])) are incorporated in the program. As the generalized- method recovers the Wood-Bosak-Zienkiewicz method (also called WBZ- method) (Wood et al.([350])), the Hilber-Hughes-Taylor method (also called HHT- method) (Hilber et al.([349])), and the Newmark family of time integration algorithms, the program allows you to take advantage of any of the these methods by specifying different input parameters, as described below.
Newmark Method
The Newmark family of time integration algorithms (Newmark([402])) is one of the most popular time integration methods as a single step algorithm. The semi-discrete equation of motion given in Equation 15–6 can be rewritten as (Hughes([166])):
(15–7) |
where:
= the nodal acceleration vector at time |
= the nodal velocity vector at time |
= the nodal displacement vector at time |
= the applied load at time |
In addition to Equation 15–7, the Newmark family of time integration algorithms requires the displacement and velocity to be updated as follows:
(15–8) |
(15–9) |
where:
= Newmark's integration parameters |
= nodal acceleration vector at time |
= nodal velocity vector at time |
= nodal displacement vector at time |
Thus, the Newmark family of time integration algorithms can be determined by the Newmark integration parameters. In the end, the Newmark integration scheme consists of the three finite difference equations presented in Equation 15–7 through Equation 15–9, and the three unknowns , , and can be numerically calculated by the three algebraic equations along with the three known quantities , , and .
By making use of the three algebraic equations given in Equation 15–7 through Equation 15–9, a single-step time integrator in terms of the unknown and the three known quantities can be written as:
(15–10) |
where:
First, the unknown is calculated using Equation 15–10. Then, the program computes the two unknowns and by using the following equations:
(15–11) |
(15–12) |
The most important factors in choosing an appropriate time integration scheme for the finite element semi-discrete equation of motion given in Equation 15–5 are accuracy, stability, and dissipation. In conditionally stable time integration algorithms, stability is affected by a chosen size of the time step; whereas in unconditionally stable time integration algorithms, a time step size can be chosen independent of stability considerations.
In the Newmark method, the amount of numerical algorithm dissipation can be controlled by one of Newmark's parameters, , as follows:
(15–13) |
With the Newmark parameters satisfying the above conditions, the Newmark family of methods may be unconditionally stable (Hughes([166])). By introducing the amplitude decay factor , the above conditions can be written:
(15–14) |
Consequently, the program provides the user with the Newmark integration procedure, which is unconditionally stable via input of the amplitude decay factor on the TINTP command. Alternatively, the and parameters may be input directly using the TINTP command.
Generalized HHT- Method
In the Newmark method, the amount of numerical dissipation can be controlled by one parameter in Equation 15–13 or in Equation 15–14. However, in low frequency modes the Newmark method fails to retain the second-order accuracy as . Note that the Newmark implicit method (constant average method; namely, and ), which is unconditionally stable and second-order accurate, has no numerical damping. If other sources of numerical damping are not introduced, the lack of numerical damping can be undesirable so that the higher frequencies of the structure can produce unacceptable levels of numerical noise (Hughes([166])).
To circumvent the drawbacks of the Newmark family of methods, the program implements the generalized HHT- method which sufficiently damps out spurious high-frequency response via introducing controllable numerical dissipation in higher frequency modes, while maintaining the second-order accuracy. It should be noted that the generalized HHT- method incorporated in the program is capable of recovering the WBZ- method (Wood et al.([350])) and the HHT- method (Hilber et al.([349])) as well as the Newmark family of time integration algorithms, depending upon the user's input on the TINTP command.
To solve for the three unknowns , , and , along with Equation 15–8 and Equation 15–9 the generalized HHT- method uses the algebraic equation:
(15–15) |
where:
Equation 15–15 give the finite difference form:
(15–16) |
where:
Analogous to the Newmark method, the generalized HHT- method calculates the unknown at time by making use of Equation 15–16. Then, the program computes the two unknowns and by using the equations given in Equation 15–11 and Equation 15–12. Since the generalized HHT- method is also an implicit time scheme, the structural stiffness matrix must be factorized to solve for at time .
As mentioned in the literature (Chung and Hulbert([348])), the generalized HHT- method is unconditionally stable and second-order accurate if the parameters meet the following conditions:
(15–17) |
where (Wood et al.([350])) and (Hilber et al.([349])). For the generalized HHT- method, you can input the four parameters on the TINTP command. By introducing the amplitude decay factor , the program also allows the user to control the amount of numerical damping if the four parameters on the TINTP command meet the following conditions:
(15–18) |
If the WBZ- method is desired, the user can control the amount of numerical damping if the four parameters on the TINTP command meet the following conditions:
(15–19) |
Finally, the program also allows a user who wants to use the generalized HHT- method to control the amount of numerical damping if the four parameters on the TINTP command meet the following conditions:
(15–20) |
It should be noted that the generalized HHT- method is second-order accurate and unconditionally stable. This method allows you to control the amount of numerical damping. The amplitude decay factor is recommended to be set as (Hughes([166])), with which any spurious participation of the higher modes can be damped out and the lower modes are not affected. A significant amount of numerical damping may be introduced by setting , but it is not recommended.
Backward Euler Method
The backward Euler method requires displacement and velocity to be updated as follows:
(15–21) |
(15–22) |
Combining equations Equation 15–21 and Equation 15–22 with the semi-discrete equation of motion given in Equation 15–7 results in a single-step time integrator in terms of the unknown {un+1} as:
(15–23) |
In nonlinear structural dynamics problems, the internal load is no longer linearly proportional to the nodal displacement, and the structural stiffness matrix is dependent on the current displacement. Therefore, Instead of Equation 15–6, any time integration scheme should be applied to the nonlinear semi-discrete equation:
(15–24) |
Equation 15–24 represents a nonlinear system of simultaneous algebraic equations; hence, any time integration operator may be used in association with the Newton-Raphson iterative algorithm. For nonlinear structural dynamics problems, both the Newmark method and the generalized HHT- method are incorporated in the program.
Newmark Method
The Newmark method assumes that at the time , the semi-discrete equation of motion given in Equation 15–24 can be rewritten as:
(15–25) |
Note that is dependent on the current displacement at time . In addition to Equation 15–25, the Newmark family of time integration algorithms requires the displacement and velocity to be updated as given in Equation 15–8 and Equation 15–9.
By introducing the residual vector , Equation 15–25 can be written as:
(15–26) |
It is important to note that the time integration operator given in either Equation 15–25 or Equation 15–26 represents a nonlinear system of simultaneous algebraic equations. Therefore, a linearized form of the time integration operator can be obtained by the Newton-Raphson method as follows:
(15–27) |
where:
= the estimate of at the kth iteration |
= the displacement increment of at the kth iteration |
Equation 15–27 gives:
(15–28) |
where:
= the tangent stiffness matrix at time |
For nonlinear structural dynamics problems, the program allows a user to input the amplitude decay factor or the Newmark integration parameters on the TINTP command.
Generalized HHT- Method
The generalized HHT- method for nonlinear structural dynamics problems assumes:
(15–29) |
where:
By introducing the residual vector , Equation 15–29 can be written as:
(15–30) |
The time integration operator given in Equation 15–29 or Equation 15–30 also represents a nonlinear system of simultaneous algebraic equations. Therefore, a linearized form of the time integration operator can be obtained by the Newton-Raphson method as follows:
(15–31) |
where:
Equation 15–31 gives:
(15–32) |
where:
Backward Euler Method
Combining the velocity and acceleration updates given in equations Equation 15–21 and Equation 15–22 with the semi-discrete equation of motion given in Equation 15–27 gives:
(15–33) |
Equation 15–33 is solved by the Newton-Raphson method. Upon convergence, the displacement is given at substep n+1, {un+1}.
Two methods of solution for the Newmark method (Equation 15–10) are available: full and mode-superposition (TRNOPT command). Each are described subsequently. Only the full solution method is available for HHT (Equation 15–15).
The full solution method (TRNOPT,FULL) solves Equation 15–10 directly and makes no additional assumptions. In a nonlinear analysis, the Newton-Raphson method (Newton-Raphson Procedure) is employed along with the Newmark assumptions. Automatic Time-Stepping discusses the procedure for the program to automatically determine the time step size required for each time step.
Inherent to the Newmark method is that the values of , , and at the start of the transient must be known. Nonzero initial conditions are input either directly (with the IC commands) or by performing a static analysis load step (or load steps) prior to the start of the transient itself. Static load steps are performed in a transient analysis by turning off the transient time integration effects (with the TIMINT,OFF command). The transient itself can then be started (by TIMINT,ON). The default with transient analysis (ANTYPE,TRANS) is for the transient to be running (TIMINT,ON); that is, to start the transient immediately. This implies . The initial conditions are outlined in the subsequent paragraphs. Cases referring to "no previous load step" mean that the first load step is transient.
- Initial Displacement
The initial displacements are:
(15–34)
where:
= vector of initial displacements = displacement vector specified by the initial conditions (IC command) = displacement vector resulting from a static analysis (TIMINT,OFF) of the previous load step - Initial Velocity
The initial velocities are:
(15–35)
where:
= vector of initial velocities = vector of velocities specified by the initial conditions (IC commands) = displacements from a static analysis (TIMINT,OFF) of the previous load step = displacement corresponding to the time point before solution. is if is the first solution of the analysis (i.e. load step 1 substep 1). = time increment between and If the previous load step was run as a static analysis (TIMINT,OFF), initial velocities are calculated using the previous two displacements and the previous time increment. Using either a single substep (NSUBST,1) or ramped loading (KBC,0) within the previous load step will result in nonzero initial velocities (assuming nonzero displacement), as shown in Figure 15.2: Effect of Number of Substeps (NSUBST) and Ramping (KBC) on Initial Velocity for TIMINT,OFF. Zero initial velocities may be obtained by using multiple substeps (NSUBST,>1) and stepped loading (KBC,1).
Figure 15.2: Effect of Number of Substeps (NSUBST) and Ramping (KBC) on Initial Velocity for TIMINT,OFF
- Initial Acceleration
The initial acceleration is simply:
(15–36)
where:
= vector of initial accelerations = vector of accelerations specified by the initial conditions (IC command). An option is available to calculate initial acceleration that is consistent with the forces applied during the start of the transient analysis. When activated by the
InitialAcc
field of the TRNOPT command, the initial acceleration is calculated as a function of applied force, inertia force, and velocity at the start time as:(15–37)
- Nodal and Reaction Load
Inertia, damping and static loads on the nodes of each element are computed.
The inertial load part of the element output is computed by:
(15–38)
where:
= vector of element inertial forces = element mass matrix = element acceleration vector The acceleration of a typical DOF is given by Equation 15–12 for time . By default, the acceleration vector is the average acceleration between time and time , since the Newmark assumptions (Equation 15–8 and Equation 15–9) assume the average acceleration represents the true acceleration. Smoothing can be suppressed using the TINTP command (AVSMOOTH option).
The damping load part of the element output is computed by:
(15–39)
where:
= vector of element damping forces = element damping matrix = element velocity vector The velocity of a typical DOF is given by Equation 15–11.
The static load is part of the element output computed in the same way as in a static analysis (Solving for Unknowns and Reactions). The nodal reaction loads are computed as the negative of the sum of all three types of loads (inertia, damping, and static) over all elements connected to a given fixed displacement node.
The mode-superposition method (TRNOPT,MSUP) uses the natural frequencies and mode shapes of a linear structure to predict the response to transient loads. This solution method imposes the following additional assumptions and restrictions:
Constant and matrices. This implies no large deflections or change of stress stiffening, as well as no plasticity, creep, or swelling.
Constant time step size.
There are no element damping matrices. However, various types of system damping are available.
Time varying imposed displacements are not allowed.
The development of the general mode-superposition procedure is described in Mode-Superposition Method. Equation 14–142 and Equation 14–143 are integrated through time for each mode by the Newmark method.
The initial value of the modal coordinates at time = 0.0 are computed by solving Equation 14–142 with and assumed to be zero.
(15–40) |
where:
= the forces applied at time = 0.0 |
The governing equation of interest is as follows:
(15–41) |
where:
= damping matrix |
= coefficient matrix |
= vector of DOF values |
= time rate of the DOF values |
= applied load vector |
In a thermal analysis, is the specific heat matrix, the conductivity matrix, the vector of nodal temperatures and the applied heat flows. Table 15.2: Nomenclature relates the nomenclature used in Derivation of Heat Flow Matrices and Derivation of Electromagnetic Matrices for thermal, magnetic and electrical analyses to Equation 15–41.
Table 15.2: Nomenclature
Thermal | temperature | heat flow |
Scalar Magnetic | scalar potential | flux |
Vector Magnetic | vector potential | current segment |
Electrical | voltage | current |
The mode-superposition procedure does not apply to first order systems.
The procedure employed for the solution of Equation 15–41 is the generalized trapezoidal rule (Hughes([166])):
(15–42) |
where:
= transient integration parameter (input
as THETA on the TINTP command) |
= nodal DOF values at time |
= time rate of the nodal DOF values at time (computed at previous time step) |
Equation 15–41 can be written at time as:
(15–43) |
Substituting from Equation 15–42 into this equation yields:
(15–44) |
The solution of Equation 15–44 employs the same solvers used for static analysis in Static Analysis. Once is obtained, is updated using Equation 15–42. In a nonlinear analysis, the Newton-Raphson method (Newton-Raphson Procedure) is employed along with the generalized trapezoidal assumption, Equation 15–42.
The transient integration parameter
(THETA
on the
TINTP command) defaults to 1.0 (backward Euler method). For
all , the system equations that follow are said to be implicit. In addition, for
the more limiting case of , the solution of these equations is said to be unconditionally stable; that
is, stability is not a factor in time step () selection. The available range of (using TINTP command) is therefore limited to
(15–45) |
which corresponds to an unconditionally stable, implicit method. is also good at damping out numerical oscillations that may result when is used.
Since influences , sudden changes in loading need to be handled carefully for values of . See the Basic Analysis Guide for more details.
The generalized-trapezoidal method requires that the values of and at the start of the transient must be known. Nonzero initial conditions are input either directly (with the IC command) (for ) or by performing a static analysis load step (or load steps) prior to the start of the transient itself. Static load steps are performed in a transient analysis by turning off the transient time integration effects (with the TIMINT,OFF command). The transient itself can then be started (TIMINT,ON). The default for transient analysis (ANTYPE,TRANS) is to start the transient immediately (TIMINT,ON). This implies . The initial conditions are outlined in Initial DOF Values.
Coupled-Physics
Some transient coupled-physics analyses involve both a first order system of transient equations and a second order system of transient equations; for example, a coupled thermal (1st order) – structural (2nd order) analysis. If the HHT time integration method is chosen in this type of analysis, the first order system of equations also follows the averaging of the forces used in the HHT time integration (described in Equation 15–15). More specifically, the governing equation for the first order system given in Equation 15–41 is modified to:
(15–46) |
where quantities with subscript follow the definitions described in Generalized HHT-α Method.
Both the Newmark and HHT methods used to solve second order equations described in Description of Structural and Other Second Order Systems are second order accurate. On the other hand, the backward Euler
method () used by default to solve the first order equations is only first order
accurate. If second order accuracy is desired for a coupled physics system such as a
piezoelectric analysis, you can request the Crank-Nicolson method () for the first order systems and constant average acceleration
( and ) by setting THETA
,
ALPHA
, and DELTA
on the
TINTP command.
Initial DOF Values
The initial DOF values for first order systems are:
(15–47) |
where:
= vector of initial DOF values |
= vector of uniform DOF values |
= DOF vector directly specified (IC command) |
= DOF vector resulting from a static analysis (TIMINT,OFF) of the previous load step available |
is set to TEMP (BFUNIF command) and/or to the temperature specified by the initial conditions (IC commands) for thermal DOFs (temperatures) and zero for other DOFs.
Nodal and Reaction Load Computation
Damping and static loads on the nodes of each element are computed.
The damping load part of the element output is computed by:
(15–48) |
where:
= vector of element damping loads |
= element damping matrix |
= element velocity vector |
The velocity of a typical DOF is given by Equation 15–42. The velocity vector is the average velocity between time and time , since the general trapezoidal rule (Equation 15–42) assumes the average velocity represents the true velocity.
The static load is part of the element output computed in the same way as in a static analysis (Solving for Unknowns and Reactions). The nodal reaction loads are computed as the negative of the sum of both types of loads (damping and static) over all elements connected to a given fixed DOF node.