Dynamic Formulation

Dynamic Timestep

The finite-difference formulation is identical to that described in Theoretical Background, except that “real” masses are used at gridpoints instead of the fictitious masses used for optimum convergence in the static solution scheme. Each tetrahedral subzone contributes one-quarter of its mass (computed from zone density and area) to each of the four associated gridpoints. The final gridpoint mass is then divided by two in the case of an eight-noded zone that contains two overlays. In finite-element terminology, FLAC3D uses lumped masses and a diagonal mass matrix.

The calculation of critical timestep is identical to that given in Theoretical Background:

(1)\[\Delta t_{\rm crit}=\min\Bigl\{ {V\over C_p \ A_{\rm max}^{f}} \Bigr\}\]

where \(C_p\) is the \(p\)-wave speed, \(V\) is the tetrahedral subzone volume, and \(A_{\rm max}^{f}\) is the maximum face area associated with the tetrahedral subzone. The min {} function is taken over all zones, and includes contributions from the structural and interface modules. A safety factor of 0.5 is used, because Equation (1) is only an estimate of the critical timestep. Hence, the timestep used for dynamic runs, \(\Delta t_d\), when no stiffness-proportional damping is used, is

(2)\[\Delta t_d=\Delta t_{\rm crit}/2\]

If stiffness-proportional damping is used (see Rayleigh Damping), the timestep must be reduced, for stability. Belytschko (1983) provides a formula for critical timestep \(\Delta t_\beta\) that includes the effect of stiffness-proportional damping:

(3)\[\Delta t_\beta=\Bigl\{{2\over\omega_{\rm max}} \Bigr\} \bigl( \sqrt{1+\lambda^2} - \lambda\bigr)\]

where \(\omega_{\rm max}\) is the highest eigenfrequency of the system, and \(\lambda\) is the fraction of critical damping at this frequency. Both \(\omega_{\rm max}\) and \(\lambda\) are estimated in FLAC3D, since an eigenvalue solution is not performed. The estimates are

(4)\[\omega_{\rm max} = {2\over\Delta t_d}\]

(5)\[\lambda={0.4\,\beta\over\Delta t_d}\]

given

(6)\[\beta = \xi_{\min}\thinspace /\thinspace \omega_{\min}\]

where \(\xi_{\min}\) and \(\omega_{\min}\) are the damping fraction and angular frequency specified for Rayleigh damping (see Rayleigh Damping). The resulting value of \(\Delta t_\beta\) is used as the dynamic timestep if stiffness-proportional damping is in operation.

Dynamic Multi-stepping

The maximum stable timestep for dynamic analysis is determined by the largest material stiffness and smallest zone in the model (see Equation (1)). Often, the stiffness and zone size can vary widely in a model (e.g., in the case of a finely zoned concrete structure located in a soft soil). A few zones will then determine the critical timestep for a dynamic analysis, even though the major portion of the model can be run at a significantly larger timestep.

A procedure known as dynamic multi-stepping is available in FLAC3D to reduce the computation time required for a dynamic calculation. In this procedure, zones and gridpoints in a model are ordered into classes of similar maximum timesteps. Each class is then run at its timestep, and information is transferred between zones at the appropriate time.

Dynamic multi-stepping uses a local timestep for each individual gridpoint and zone. At the start of an analysis, the grid is scanned and the local stable timestep for each gridpoint, \(\Delta t_{gp}\), is determined and stored. The value of \(\Delta t_{gp}\) depends on size, stiffness and mass of the neighboring subzones (as shown in Equation (1)), attached structural elements and interfaces. The global timestep, \(\Delta t_{G}\), is determined as the minimum of all \(\Delta t_{gp}\), as in the standard formulation.

Integer multipliers, \(M_{gp}\), to the global timestep are then determined for each gridpoint according to the algorithm illustrated by the flow chart in Figure 1. This algorithm ensures that multipliers are powers of 2. In the current implementation, \(M_{gp}\) is set to 1 for nodes that are assigned a null material model, connected to structural elements, attached to other gridpoints, or part of a quiet boundary. All zones are then scanned, and an integer multiplier, \(M_{z}\), is calculated for each zone as the minimum of the multipliers for the four surrounding gridpoints.

Figure 1: Flow chart for determination of gridpoint multiplier, \(M_{gp}\).

Calculations for a zone (i.e., derivation of new stresses from surrounding gridpoint velocities; accumulation of gridpoint force sums from stress components) are only performed every \(M_{z}\) timesteps. In all expressions involving a timestep, the global timestep is replaced by \(\Delta t_{G}\) \(M_{z}\).

Calculations for a gridpoint (i.e., derivation of new velocities and displacements from gridpoint force sums) are only performed every \(M_{gp}\) timesteps; otherwise, the force sums are reset to zero, which is normally done after every motion calculation. In all expressions involving a timestep, the global timestep is replaced by \(\Delta t_{G}\) \(M_{gp}\).

The effect of the prescriptions described above is to skip calculation of selected gridpoints and zones, thereby speeding up the overall calculation. The use of gridpoint and zone multipliers (\(M_{gp}\) and \(M_{z}\), respectively) ensures the following characteristics:

1. The force sum at each gridpoint is composed of component forces from each connected zone that exist at the same point in time. The simultaneous nature of the component forces is guaranteed by the fact that multipliers are powers of two. Arbitrary integral multipliers would not have this characteristic.
2. Velocities seen by a zone (at the four surrounding gridpoints) are not updated between zone updates. This is guaranteed by the fact that the zone multiplier is the minimum of the surrounding gridpoint multipliers. Since stress increments are derived from strain and displacement increments, the displacement contribution of a gridpoint is felt by a zone at each update, even though the gridpoint is updated less frequently than the zone. In essence, the total displacement increment of the gridpoint is divided into \(M_{gp}/M_{z}\) equal parts.

This scheme is accurate for dynamic simulations that represent waves with frequencies well below the natural frequencies of individual elements. The condition is usually guaranteed by the wavelength criterion, \(\Delta l \le {\lambda / 10}\), which is described in Accurate Wave Propagation. For higher frequencies, it is believed that inaccuracies arise from the fact that velocities used in computing strain increments are not defined (in time) at the center of the time interval, \(\Delta t\), for the case of a zone multiplier being unequal to the gridpoint multiplier. This represents a departure from the second-order accuracy of the central difference scheme used in FLAC3D. However, it is always possible to assess the accuracy of the scheme for any part of the simulation by running a short period of the simulation with and without dynamic multi-stepping; the results may be directly compared.

Dynamic multi-stepping is invoked with the command zone dynamic multi-step on. When cycling begins, the division of zones into the different multiplier classes will be printed to the screen. The effect of dynamic multi-stepping on calculation speed is model-dependent (i.e., the more zones that have a high multiplier, the greater the increase in speed).

The following example is an illustration of the effects of dynamic multi-stepping:

• Shear Wave On Stiff Wall In Soft Soil: Dynamic Multi-Stepping

The multi-stepping multipliers can be plotted as a label plot on zones; the gridpoint multipliers can be contoured.

Dynamic multi-stepping can be used with structural elements.

For additional information and example applications of dynamic multi-stepping, see Unterberger, Cundall, and Zettler (1997). The application of dynamic multi-stepping in numerical predictions of vibrations caused by rail traffic in tunnels is presented in Unterberger, Hochgatterer, and Poisel (1997) and Daller, Unterberger, and Hochgatterer (1996).