The automatic time-stepping (or automatic loading) method adjusts the time-step size and/or the applied loads automatically in response to the current state of the analysis. This method (accessed via AUTOTS,ON) can be applied to structural, thermal, electric, and magnetic analyses performed in the time domain (TIME), and includes static (or steady-state) (ANTYPE,STATIC) and dynamic (or transient) (ANTYPE,TRANS) situations.
An important point to be made here is that automatic loading always works through the adjustment of the time-step size; and that the loads that are applied are automatically adjusted if ramped boundary conditions are activated (using KBC,0). In other words the time-step size is always subjected to possible adjustment when automatic loading is engaged. Applied loads and boundary conditions, however, will vary according to how they are applied and whether the boundary conditions are stepped or ramped. That is why this method may also be thought of as automatic loading.
The automatic time-stepping algorithm has two key features:
The capability to estimate the next time-step size. The esitmation is based on current and past analysis conditions, and the algorithm makes proper load adjustments as necessary. That is, given conditions at the current time, tn, and the previous time increment, Δtn, the primary aim is to determine the next time increment, Δtn+1. Since the determination of Δtn+1 is largely predictive, this part of the automatic time-stepping algorithm is referred to as the time-step prediction component.
The time-step bisection component. Its purpose is to decide whether or not to reduce the present time-step size, Δtn, and redo the substep with a smaller step size. For example, working from the last converged solution at time point tn-1, the present solution begins with a predicted time step, Δtn. Equilibrium iterations are performed; and if proper convergence is either not achieved or not anticipated, this time step is reduced to Δtn/2 (that is, it is bisected), and the analysis begins again from time tn-1. Multiple bisections can occur per substep for various reasons (discussed later).
At a given converged solution at time tn, and with the previous time increment Δtn, the goal is to predict the appropriate time-step size to use as the next substep. The predicted time-step size is derived from the results of several unrelated calculations and is most easily expressed as the minimization statement:
(14–77) |
where:
Δteq = time increment limited by the number of equilibrium iterations needed for convergence at the last converged time point. The more iterations required for convergence, the smaller the predicted time step. This is a general measure of all active nonlinearities. Increasing the maximum number of equilibrium iterations (NEQIT) will tend to promote larger time-step sizes. |
Δt1 = time increment limited by the response eigenvalue computation for first-order systems (such as thermal transients) (TINTP). |
Δt2 = time increment limited by the response frequency computation for second-order systems (such as structural dynamics). The aim is to maintain 20 points per cycle (described below). When the middle step criterion is used, this criterion can be disabled. |
Δtg = time increment representing the time point at which a gap or a nonlinear (multi-status) element will change abruptly from one condition to another (status change). KEYOPT(7) enables further control for the contact elements. |
Δtc = time increment based on the allowable creep-strain increment. |
Δtp = time increment based on the allowable plastic-strain increment. The limit is set at 5 percent per time step. |
Δtm = time increment limited by the middle step residual tolerance (described below) for second-order systems (such as structural dynamics) (MIDTOL). When enabled, the Δt2 criterion can be disabled. |
The program calculates several trial step sizes and selects the minimum step size for the next time step. The predicted value is further restricted to a range of values expressed by:
(14–78) |
and
(14–79) |
where:
F = increase/decrease factor. F = 2, if static analysis; F = 3, if dynamic (see the ANTYPE and TIMINT commands) |
Δtmax = maximum time-step size (DTMAX on the DELTIM command or the equivalent quantity calculated from the NSUBST command) |
Δtmin = minimum time-step size (DTMIN on the DELTIM command or the equivalent quantity calculated from the NSUBST command) |
In other words, the current time-step is increased or decreased by a factor of 2 (or 3 if dynamic) at most, and cannot be less than Δtmin or greater than Δtmax.
It is also possible to control the time-step cutback during a nonlinear solution (CUTCONTROL), and to set a decision parameter for automatically increasing the time-step interval (OPNCONTROL).
When bisection occurs, the current substep solution (Δtn) is removed, and the time-step size is reduced by 50 percent. If applied loads are ramped (KBC,0), then the current load increment is also reduced by the same amount. One or more bisections can take place for several reasons, namely:
The number of equilibrium iterations used for this substep exceeds the number allowed (NEQIT command).
It appears likely that all equilibrium iterations will be used.
A negative pivot message was encountered in the solution, suggesting instability.
The largest calculated displacement DOF exceeds the limit (DLIM on the NCNV command).
An illegal element distortion is detected (for example, negative radius in an axisymmetric analysis).
For transient structural dynamics, when the middle step residual is greater than the given tolerance. This check is done only when the middle step residual check is enabled by the MIDTOL command.
More than one bisection may be performed per substep. However, bisection of the time-step size is limited by the minimum size (defined by DTMIN input on the DELTIM command or the equivalent NSUBST input).
The response eigenvalue is used in the computation of Δt1 and is defined as:
(14–80) |
where:
λr = response eigenvalue (item RESEIG for POST26 SOLU command and *GET command) |
{Δu} = substep solution vector (tn-1 to tn) |
[KT] = the Dirichlet matrix. In a heat transfer or an electrical conduction analysis this matrix is referred to as the conductivity matrix; in magnetics this is called the magnetic "stiffness". The superscript T denotes the use of a tangent matrix in nonlinear situations |
[C] = the damping matrix. In heat transfer this is called the specific heat matrix. |
The product of the response eigenvalue and the previous time step (Δtn) has been employed by Hughes([146]) for the evaluation of 1st order explicit/implicit systems. In Hughes([146]) the quantity Δtnλ is referred to as the "oscillation limit", where λ is the maximum eigenvalue. For unconditionally stable systems, the primary restriction on time-step size is that the inequality Δtnλ >> 1 should be avoided. Hence it is very conservative to propose that Δtnλ = 1.
Since the time integration used employs the trapezoidal rule (Equation 15–42), all analyses of 1st order systems are unconditionally stable. The response eigenvalue supplied by means of Equation 14–80 represents the dominate eigenvalue and not the maximum; and the time-step restriction above is restated as:
(14–81) |
This equation expresses the primary aim of automatic time-stepping for 1st order transient analyses. The quantity Δtnλr appears as the oscillation limit output during automatic loading. The default is f = 1/2, and can be changed (using OSLM and TOL on the TINTP command). The quantity Δt1 is approximated as:
(14–82) |
If the thermal solution reaches steady-state, the time step is allowed to increase by a factor of 2. Steady-state is determined if the maximum temperature change is less than 3 degrees for 3 time steps in a row. Use the OPNCONTROL command to modify these defaults, if desired.
For the Quasi option (THOPT command), the time step is also controlled by the property change tolerance such that the chosen time step keeps the property change below the tolerance. See the THOPT command for more details.
The response frequency is used in the computation of Δt2 and is defined as (Bergan([106])):
(14–83) |
where:
fr = response frequency (item RESFRQ for POST26 SOLU command and *GET command) |
{Δu} = substep solution vector (tn-1 to tn) |
[KT] = tangent stiffness matrix |
[M] = mass matrix |
This equation is a nonlinear form of Rayleigh's quotient. The related response period is:
(14–84) |
Using Tr the time increment limited by the response frequency is:
(14–85) |
When the middle step criterion is used, this criterion can be turned off.
The time-step size may be increased or decreased by comparing the value of the creep ratio Cmax (Rate-Dependent Plasticity (Including Creep and Viscoplasticity)) to the creep criterion Ccr. Ccr is equal to .10 unless it is redefined (using the CRPLIM command). The time-step estimate is computed as:
(14–86) |
Δtc is used in Equation 14–77 only if it differs from Δtn by more than 10%.
The time-step size is increased or decreased by comparing the value of the effective plastic strain increment (Equation 4–26) to 0.05 (5%). The time-step estimate is computed as:
(14–87) |
Δtp is used in Equation 14–77 only if it differs from Δtn by more than 10%.
The midstep residual is used in the computation of Δtm. The midstep residual for the determination of the time step is based on the following consideration. The solution of the structural dynamic analysis is carried out at the discrete time points, and the solution at the intermediate time remains unknown. However, if the time step is small enough, the solution at the intermediate time should be close enough to an interpolation between the beginning and end of the time step. If so, the unbalanced residual from the interpolation should be small. On the other hand, if the time step is large, the interpolation will be very different from the true solution, which will lead to an unbalanced residual that is too large. The time step is chosen to satisfy the criterion set by the user (e.g. MIDTOL command).
Refer to the discussion in Newton-Raphson Procedure. The residual force at any time between the time step n and n+1 can be written as:
(14–88) |
where:
ν = intermediate state between the time step n and n+1 (0 < ν < 1) |
{R} = residual force vector |
= vector of the applied load at n + ν |
= vector of the restoring load corresponding to the element internal load at n + ν, which depends on the intermediate state of displacement at n + ν, and also the velocity and acceleration at n + ν. This intermediate state is approximately calculated based on the Newmark assumption. |
A measure of the magnitude of {R} is established in a manner similar to the convergence check at the end of the time step. (See Convergence). After the solution has converged at the end of the time step (n+1), the midstep residual force is compared to the reference value:
(14–89) |
where:
||{R}|| = magnitude (vector norm) of residual force vector |
Rref = reference force (see Convergence) |
The convergence criterion for the midstep residual is defined by the value of
τb (input as TOLERB
on
MIDTOL command):
If τb > 0, the value is used as a tolerance. If τb = 0 is specified or τb is not specified, then a default positive value is used as a tolerance. The midstep residual is assumed to have converged if its value is within the desired tolerance (ε τb ). Depending on how well the convergence criterion is satisfied the time-step size for the next increment is increased or kept unchanged.
If the midstep residual hasn’t converged (ε > τb), the time step is repeated with a smaller increment:
(14–90) |
where:
Δtn = old time-step size |
τb = midstep residual
tolerance(TOLERB on MIDTOL
command) |
If τb < 0, the value is used as a reference force (reference moment is computed from reference force value) for midstep convergence check. A procedure similar to the one described above is followed with modified definition of time-step size:
(14–91) |