Nonlinear Transient Analysis

Nonlinear Small Displacement Transient and Nonlinear Large Displacement Transient (LGDISP) Analysis are currently available to solve nonlinear problems which include transient effects.

The predominant difference between Nonlinear Static Analysis and Nonlinear Transient Analysis is the inclusion of inertia and momentum terms in the solution of the equation of motion. This nonlinear transient solution sequence typically supports all nonlinear features supported by Nonlinear Static Analysis, including geometric nonlinearity, material nonlinearity, and contact. Subcase continuation between Nonlinear Transient Analysis and Nonlinear Static Analysis is also supported.

Nonlinearity Sources

Geometric Nonlinearity

In analyses involving geometric nonlinearity, changes in geometry as the structure deforms are considered in formulating the constitutive and equilibrium equations. Many engineering applications require the use of large deformation analysis based on geometric nonlinearity. For example, applications such as metal forming, tire analysis, and medical device analysis. Small deformation analysis based on geometric nonlinearity is required for some applications, like analysis involving cables, arches and shells. Such applications involve small deformation, except finite displacement or rotation.

Material Nonlinearity

Material nonlinearity involves the nonlinear behavior of a material based on current deformation, deformation history, rate of deformation, temperature, pressure, and so on.

Constraint and Contact Nonlinearity

Constraint nonlinearity in a system can occur if kinematic constraints are present in the model. The kinematic degrees-of-freedom of a model can be constrained by imposing restrictions on its movement. In OptiStruct, constraints are enforced with Lagrange multipliers.

In the case of contact, the constraint condition is based on inequalities and such a constraint generally does not allow penetration between any two bodies in contact.

Follower Load

Applied loads can depend upon the deformation of the structure when large deformations are involved. Geometrically, the applied loads (Forces or Pressure) can deviate from their initial direction based on how the model deforms at the location of application of load. In OptiStruct, if the applied load is treated as follower load, the orientation and/or the integrated magnitude of the load will be updated with changing geometry throughout the analysis.

Nonlinear Transient Solution Method

Nonlinear Transient solutions are done by marching through the solution in the time domain. The time-discretized equation is solved using Newton's method (similar to NLSTAT). Two time stepping schemes are provided for nonlinear transient analysis. The TMTD field on the TSTEP entry can be used to switch between the two integration schemes. The Generalized Alpha method (TMTD=1, default) and the Backward Euler method (TMTD=2).

Generalized Alpha (α) Method

In the Generalized Alpha method, the equilibrium equation takes the following form:(1)

$f_{e x t}^{t + α h} = M ((1 - α_{m}) a^{t + α h} + α_{m} a^{t}) + C v^{t + α h} + f_{i n t}^{t + α h}$

Where,

$M$: Mass
$C$: Viscous damping
$K$: Stiffness matrix
$f$: Total force

Subscripts, $e x t$ and $i n t$ indicating external and internal force, respectively. Superscripts, $t$ and $t + h$ indicate the time when the variable is calculated, with $t$ being the previous time and $t + h$ being the current time when the displacement increment is being solved for. The Lagrangian coordinate is denoted by $x$ , and the material time derivatives of $x$ are denoted by $v$ (velocity) and $a$ (acceleration).

The velocity and acceleration at current time are computed with Newmark approximation by involving parameters

$γ$ and

$β$ .(2)

$v^{t + h} = v^{t} + h ((1 - γ) a^{t} + γ a^{t + h})$

(3)

$u^{t + h} = u^{t} + h v^{t} + \frac{1}{2} h^{2} ((1 - 2 β) a^{t} + 2 β a^{t + h})$

The superscript

$t + α h$ indicates, for a general quantity,

$z$ , that is using parameter

$α$ :(4)

$z^{t + α h} = (1 + α) z^{t + h} - α z^{t}$

As the only special case, the acceleration is interpolated using parameter

$α_{m}$ .(5)

$a^{t + α h} = (1 + α_{m}) a^{t + h} - α_{m} a^{t}$

From the above equations, it can be inferred that the Generalized Alpha method is controlled by four non-dimensional parameters (

$α, β, γ,$ and

$α_{m}$ ). With

$α_{m} = 0.0$ or

$α = α_{m} = 0.0$ , the method degenerates into the HHT-

$α$ method or Newmark-

$β$ method, respectively. The parameters can be selected as:(6)

$[- \frac{1}{3} < α \leq 0]; [β \geq \frac{1}{4} - \frac{1}{2} (α_{m} + α)]; [- 1 \leq α_{m} < \frac{1}{2}]$

Given

$α_{m}$ and

$α$ , the default values of

$γ$ and

$β$ which ensure the method is unconditionally stable (for linear problems), 2^nd order accurate and with maximized high frequency dissipation:(7)

$[γ = \frac{1}{2} - (α_{m} + α)]; [β = \frac{1}{4} {(1 - α_{m} - α)}^{2}]$

By default, $α_{m} = 0.0$ and $α = - 0.05$ . That is, the default scheme is HHT- $α$ method.

The Generalized Alpha method is solved using Newton's method. In each iteration,

$Δ u^{t + h}$ , the increment of displacement

$Δ u^{t + h} = u^{t + h} - u^{h}$ , is solved for. The subscript denotes the iteration count.(8)

$K^{*} Δ u_{j + 1}^{t + h} = f_{e x t}^{t + h} - f_{i n t}^{t + α h} - M a_{j}^{t - α_{m} h} - C v_{j}^{t + h}$

(9)

$K^{*} = \frac{1 - α_{m}}{β h^{2}} M + (1 + α) \frac{γ}{h β} C + (1 + α) K_{j}$

(10)

$a_{j}^{t - α_{m} h} = α_{m} a^{t} + (1 - α_{m}) a_{j}^{t + h} = α_{m} a^{t} + (1 - α_{m}) (\frac{1}{β h^{2}} Δ u_{j}^{t + h} - \frac{1}{β h} v^{t} - \frac{1 - 2 β}{2 β} a^{t})$

(11)

$v_{j}^{t + h} = v^{t} + h ((1 - γ) a^{t} + γ a_{j}^{t + h})$

Backward Euler Method

In the Backward Euler method, the equilibrium equation takes the following form:(12)

$f_{e x t}^{t + α h} = M a^{t + h} + C v^{t + h} + f_{i n t}^{t + α h}$

(13)

$v^{t + h} = \frac{u^{t + h} - u^{t}}{h}$

(14)

$a^{t + h} = \frac{v^{t + h} - v^{t}}{h}$

Where,

$M$: Mass
$C$: Viscous damping
$K$: Stiffness matrix
$f$: Total force

The Backward Euler method is also solved using Newton's method. In each iteration,

$Δ u^{t + h}$ , the increment of displacement

$Δ u^{t + h} = u^{t + h} - u^{h}$ , is solved for. The subscript denotes the iteration count.(15)

$K^{*} Δ u_{j + 1}^{t + h} = f_{e x t}^{t + h} - f_{i n t}^{t + h} - M a_{j}^{t + h} - C v_{j}^{t + h}$

(16)

$K^{*} = \frac{1}{h^{2}} M + \frac{1}{h} C + K_{j}$

(17)

$a_{j}^{t + h} = \frac{1}{h^{2}} Δ u_{j}^{t + h} - \frac{1}{h} v^{t}$

(18)

$v_{j}^{t + h} = \frac{1}{h} Δ u_{j}^{t + h}$

The only parameter for Backward Euler method is the time step $h$ .

The Backward Euler method does not require the input of TCi fields on the TSTEP entry. The Alpha and Beta fields introduce subcase-dependent Rayleigh damping, so that the viscous damping matrix $C$ in a particular subcase is $C = α M + β K$ .

Typically, the Generalized Alpha method should be used for most Nonlinear Transient Analysis. In this method, numerical damping can be adjusted through the parameters $α$ , $β$ , $γ$ and $α_{m}$ . In particular, non-zero $α$ and $α_{m}$ introduce damping for high-frequency response components. However, Backward Euler method can be used for Static Analysis such as post-buckling solutions, since this method is dissipative and therefore highly stable.

Automatic Time Stepping

This Nonlinear Transient solution provides automatic time stepping based on Local Truncation Error (LTE).

The HHT-

$α$ method is used as an example to derive the method, and the other methods follow similar steps. Using HHT-

$α$ method, by Taylor's expansion, the LTE of displacement

$u$ from time

$t$ to

$t + h$ is:(19)

$\begin{array}{l} e = \frac{1}{6} h^{3} \overset{⃛}{u} (τ) \\ t < τ < t + h \end{array}$

Using acceleration,

$e$ can be approximated as:(20)

$e = \frac{1}{6} h^{2} |Δ a|$

For coupled multiple DOF system, the above absolute value is replaced by certain norm of

$Δ a$ . Here the mass-weighted norm is used, so that:(21)

$e = \frac{1}{6} h^{2} {(Δ a^{T} M Δ a)}^{\frac{1}{2}}$

Error estimation in the previous equation requires normalization, as it depends on certain reference displacement, for example, initial displacement, as well as cut-off frequency. The normalized error estimation is:(22)

$\tilde{e} = \frac{e}{K_{Ω} ‖ u_{R} ‖}$

Here,

$K_{Ω}$ is the time-averaged Local Truncation Error (LTE) measure of a linear undamped single DOF oscillator with unit initial displacement and oscillating at approximately 10 cycles per second or with cut-off frequency

$Ω = ω h = 0.6$ (where,

$ω$ is the natural frequency of the single DOF system).

$‖ u_{R} ‖$ is the maximum value of the displacement norm over all the previous time steps.(23)

$‖ u_{R} ‖ = \max_{1 \leq j \leq n} (‖ u_{j} ‖)$

Where, $u_{j}$ is the displacement at time step, $j$ .

The norm is defined as:(24)

$‖ u_{j} ‖ = \sqrt{{(u_{j})}^{T} M (u_{j})}$

The time steps are automatically adjusted based on the following conditions (TOL is the user-defined tolerance set on the TSTEP Bulk Data Entry), if:

$\tilde{e}$ > TOL - Reject current step, cutback and redo the current step.
TOL > $\tilde{e}$ > 0.5 * TOL - Accept current step, cutback the next step.
0.5 * TOL > $\tilde{e}$ > 1/16 * TOL - No changes.
1/16 * TOL > $\tilde{e}$ - The next time step is doubled.

The MREF continuation line on TSTEP entry can be used to control automatic time stepping, so that the time step $h$ is adjusted according to the LTE of the current step. As shown above, when error is "large" when compared to the tolerance (TOL), $h$ will be reduced by half and the current step is re-calculated. The maximum number of such operations within each step is controlled by the TN1 field. On the other hand, when $h$ is "small" compared to the tolerance (TOL), $h$ is requested to be increased, but only after TN2 contiguous steps with such a request, will $h$ be actually increased.

Damping

This Nonlinear Transient solution sequence currently only provides support for Rayleigh Damping. The global Rayleigh damping parameters can be input using PARAM, ALPHA1 and PARAM, ALPHA2. Subcase-dependent Rayleigh damping parameters can be input using the Alpha and Beta fields on the TSTEP Bulk Data Entry. Material-dependent Rayleigh damping parameters can be input via the RAYL continuation line on the corresponding material entry (MAT1, MAT2, MAT8, MAT9, and MAT9OR Bulk Data Entries).

Baseline Correction

Baseline correction can be applied to correct the acceleration record to avoid displacement shift.

When acceleration is applied as input via SPCD, the integration of displacement results may diverge at the last point of simulation. Baseline correction can be used to correct the acceleration input. The BASELIN Subcase and BASELIN Bulk Data Entry pair can be used to apply baseline correction.

There are two methods by which Baseline correction can be applied. For both cases, the acceleration load is corrected by fitting the load to a polynomial with coefficients $C_{0}$ , $C_{1}$ , and $C_{2}$ .

Acceleration Method (METHOD=ACCE on BASELIN Bulk Data)

The baseline correction to the acceleration load is calculated as: (25)

$a_{b} (t) = C_{0} + C_{1} t + C_{2} t^{2}$

The corrected acceleration load is:(26)

$a_{c} (t) = a (t) - a_{b} (t)$

The coefficients of the polynomial are evaluated by minimizing the square of the corrected acceleration load.(27)

$S (C_{1}, C_{2}, C_{3}) = \min [\sum_{i = 1}^{n} {(a_{c} (t_{i}))}^{2}]$

Where,

$C_{0}$ , $C_{1}$ , and $C_{2}$: Coefficients of the polynomial, which are unknowns to be calculated from the above minimization functions.
$S$: Function of polynomial coefficients.
$a (t)$: Input acceleration load.
$a_{b} (t)$: Baseline correction of the acceleration.
$a_{c} (t)$: Corrected acceleration after baseline correction.
$t$: Time associated with each acceleration data point in the record.
$n$: Total number of acceleration data points in the record (points on TABLED# entry referenced on TLOAD1 entry).

Velocity Method (METHOD=VELO on BASELIN Bulk Data)

The baseline correction to the acceleration load is calculated as: (28)

$a_{b} (t) = C_{0} + C_{1} (\frac{t - T_{1}}{T_{2} - T_{1}}) + C_{2} {(\frac{t - T_{1}}{T_{2} - T_{1}})}^{2}$

The corrected acceleration load is:(29)

$a_{c} (t) = a (t) - a_{b} (t)$

The corrected velocity is:(30)

$v_{c} (t) = v (T_{1}) + \int_{T_{1}}^{T_{2}} a_{c} (t) d t$

The coefficients of the polynomial are evaluated by minimizing the integration of square of the corrected velocity load. (31)

$S (C_{1}, C_{2}, C_{3}) = \min [\int_{T_{1}}^{T_{2}} {(v_{c} (t))}^{2} d t]$

Where,

$C_{0}$ , $C_{1}$ , and $C_{2}$: Coefficients of the polynomial, which are unknowns to be calculated from the above minimization functions.
$S$: Function of polynomial coefficients.
$a_{b} (t)$: Baseline correction of the acceleration.
$a_{c} (t)$: Corrected acceleration after baseline correction.
$v (t)$: Velocity load.
$v_{c} (t)$: Corrected velocity after baseline correction.
$t$: Time associated with each acceleration data point in the record.
$T_{1}$ and $T_{2}$: Initial and final total time.; For both continuation subcases or independent subcases, the baseline correction calculation is conducted for each subcase, and separate $C_{0}$ , $C_{1}$ , and $C_{2}$ coefficients are calculated for each subcase for baseline correction. In the case of continuation subcases, the values of $T_{1}$ and $T_{2}$ used in the calculation are still the total time.

Problem Setup

Input

For Nonlinear Transient Analysis, the loads and boundary conditions are defined in the Bulk Data Entry section of the input deck. The static loads and boundary conditions are defined in the Bulk Data Entry section of the input deck. They should be referenced in the Subcase Information Entry section using a SPC and DLOAD statement in a subcase, respectively.

For Nonlinear One-step Transient Thermal Stress Analysis (OSTTS), if the Nonlinear Transient Analysis subcase has a temperature load defined (via TEMPT), then a TEMP Subcase Information Entry referencing the TEMPT Bulk Data Entry can be directly specified in the nonlinear transient subcase.

Large Displacement Analysis can be activated via PARAM, LGDISP, 1, NLPARM(LGDISP), or NLCTRL(LGDISP) Subcase Information Entries.

To indicate that a nonlinear transient solution is required for any subcase, the combination of NLPARM (or NLCTRL) and TSTEP entries should be present (an alternative definition is to just use TSTEPNL instead of both NLPARM (or NLCTRL) and TSTEP). The NLPARM (or NLCTRL) Subcase Information Entry, in turn, points to the NLPARM (or NLCTRL) Bulk Data Entry, and the TSTEP entry points to the TSTEP Bulk Data Entry. If TSTEPNL is used in the subcase, then it points to the TSTEPNL Bulk Data Entry.

Using both TSTEPNL and NLPARM, or TSTEPNL and NLCTRL in a subcase is not supported (refer to the TSTEPNL Bulk Data Entry for more information). If NLPARM or NLCTRL is present in a Nonlinear Transient subcase, then a TSTEP entry is mandatory. However, if TSTEPNL entry is present instead of NLPARM or NLCTRL, then the TSTEP entry is not compulsory.

The NLOUT Bulk Data Entry and NLOUT Subcase Information Entry can be used to control incremental output. The NLADAPT Bulk Data Entry and Subcase Information Entry can be used to define parameters for time-stepping and convergence criteria.

The DLOAD Subcase Entry can reference corresponding DLOAD or TLOAD1/TLOAD2 Bulk Data Entries. The TLOAD1/TLOAD2 entries define time-dependent loading wherein the EXCITEID field can point to applied loads (DAREA, FORCEx, PLOADx, etc.), enforced motion (SPC/SPCD) or temperature loads (TEMP/TEMPD).

▸Example (NLPARM)

SUBCASE 10
   ANALYSIS = DTRAN
   SPC = 1
   DLOAD = 2
   NLPARM = 99
   TSTEP = 2
   NLOUT = 23
   IC = 12
.
.
BEGIN BULK
PARAM,LGDISP,1
NLPARM  99
TSTEP, 2

▸Example (NLCTRL)

SUBCASE 10
   ANALYSIS=DTRAN
   SPC = 1
   DLOAD = 2
   NLCTRL = 99
   TSTEP = 2
   NLOUT = 23
   IC = 12
.
.
BEGIN BULK
PARAM,LGDISP,1
NLCTRL  99 
TSTEP, 2

Subcase Continuation

Subcase continuation is supported between Nonlinear Transient subcases in the same model. It is also supported from a NLSTAT subcase to a Nonlinear Transient subcase.

With subcase continuation, subcase 2 follows subcase 1. When both are of Nonlinear Transient type, the initial condition (IC entry, if present) in subcase 2 is ignored. Instead, subcase 2 obtains its initial condition from the last time step of subcase 1.

On the other hand, if subcase 1 is of NLSTAT type, IC entry is allowed for subcase 2. In this situation you must also be cautious in defining the loading of subcase 1. If LOAD is used to define loading for subcase 1, its magnitude will be gradually ramped down to zero at the beginning of subcase 2. If DLOAD is used in subcase 1, the load will not continue into subcase 2. Instead it becomes zero immediately when entering subcase 2. In other words, for NLSTAT(LOAD) + Nonlinear Transient, the load in the 2nd subcase is the combination of DLOAD of Nonlinear Transient and the (ramping down) LOAD of NLSTAT. Whereas, for NLSTAT(DLOAD) + Nonlinear Transient, the load in the 2nd subcase is only from DLOAD of Nonlinear Transient.

One typical use case for subcase continuation is the post-buckling analysis. One can load the structure to a point close to buckling using a NLSTAT Subcase, followed by a Nonlinear Transient subcase (using the Backward Euler method to improve stability) for the rest of buckling process.

Subcase continuation between Small and Large Displacement Nonlinear Transient Analyses (or vice-versa) is not allowed.

Output

The typical output entries (DISPLACEMENT, VELOCITY, and ACCELERATION) can be used to request corresponding output for Nonlinear Transient Analysis. The NLOUT Subcase and Bulk Data Entries can be used to request intermediate results.

Energy output is available by default and additional control is available via NLENRG Bulk/Subcase pair.

If automatic time-stepping is not active, then the typical convergence results are output for each time step similar to NLSTAT Subcases. However, if automatic time-stepping is activated, then, in addition to the convergence results, an extra line is added to each time step.

The fields show the current time step number (itstp), current time (time), norm of acceleration error (norm(da)), the norm of the reference displacement (u_ref), the acceleration error (err_da) and adjustment request (adjustment), respectively.

The adjustment request field shows the following syntax:

Adjustment Printing (.out file): Information
Reduce Next: Indicates a reduction of next time step for $h$
Enlarge Next: Indicates an increase of next time step for $h$
No Change: Indicates no change in $h$ for current and subsequent step.
Cutback: Indicates reduction of current step. The current step is then redone.

Note:

The same limitations that apply to Nonlinear Static Analysis also apply to Nonlinear Transient Analysis.
The NLADAPT entry is recommended to be used to control the maximum and minimum time-step in Nonlinear Transient Analysis.
Linear and large rotation DMIG (Direct Matrix Input at Grid points) are supported in Nonlinear Transient Analysis.

Plot Convergence During Analysis

The convergence plot can be obtained for all nonlinear iterations across all subcases for nonlinear static and nonlinear transient analysis, using the Altair Compute Console (ACC).

The progress can also be monitored for all the optimization iterations.

For more details, refer to Plot Convergence During Analysis in the Runtime Monitoring section.