Nodal Time Step Subcycling

Using a uniform time step in non-uniform meshes is a very costly and inefficient. The extensive research on the topic of time step optimization has led to the development of various multi-time step subcycling algorithms. The initial work of Belytschko and Mullen in 1976 led to an ``implicit-explicit'' method for second order equations using nodal partitioning (Neal and Belytschko, 1998). Further research, by Hughes and Lui (1978), led to a ``implicit-explicit'' subcycling method using element partitioning. In our research we have adopted an explicit subcycling method for second order equations presented by Smolinski (1989). This algorithm allows the use of different time steps in different regions of the same problem. As a result, we are not constrained by the minority elements, rather each region is given a time step that maximizes its efficiency while still satisfying the local Courant condition. This has the added benefit of minimizing the number of calculations necessary in the larger time step regions over the same time period. Multi-time step methods partition a mesh with different time steps using either nodal or element partitioning. The algorithm presented by Smolinski employs nodal partitioning. For nodal partitioning, time steps are distributed to nodes or nodal groups, while in element partitioning, the elements themselves are given different time steps.

Although the algorithm presented by Smolinski is capable of supporting multiple time step regions, we have decided, for simplicity reason, to limit the number of regions to only two - designated as regions $A$ and $B$. It should be noted that the nodes of the two regions are not restricted to be grouped together, instead they can be dispersed over the entire mesh as is shown in Figure 2.5.

Figure 2.5: Subcycling region distribution.

For our analysis, we designate region $A$ as the critical region, having the smaller time step equal to $1 \Delta t$. This region encompasses the smaller volumetric elements as well as any cohesive elements in the system. Region $B$ is given a time step of $m \Delta t$, where the subcycling parameter, $m$, is thus the ratio of time steps from region $B$ to region $A$. If $m = 1$, all regions are given an equal time step and no subcycling is performed. When subcycling is used, this time step ratio must always be a positive even number.

Once the nodal partitioning has been performed for the initial mesh, the subcycling algorithm is directly applied to the central difference time stepping loop. As the solution is stepped from time $t$ to time $t+m \Delta t$, the nodes of region $A$ are updated $m$ times using the standard discretization equations. Region $B$, on the other hand, is the subcycled region over the same time period, to which an approximate method is applied. Region $B$ is not explicitly updated, rather the nodal accelerations of region $B$ are toggled after each time step. Since the parameter $m$ is even and since the nodal velocity update is of the form ${\bf v_{n+1}} = {\bf v_{n}} + f({\bf a}_n + {\bf a}_{n+1})$, alternating the sign of the acceleration has for effect to keep the velocity constant in region $B$. This toggling of nodal accelerations in region $B$ foregoes the calculation of the internal and cohesive forces, thereby saving a significant portion of time that would normally be dedicated to these calculations. The equations used for direct calculations in region $A$ and approximations in region $B$ are thus given by

$${\bf d}^A_{n+1} = {\bf d}^A_n + \Delta t {\bf v}^A_n + \frac{1}{2} \Delta t^2 {\bf a}^A_n,\vspace*{0.25cm}$$ (2.16)

$${\bf d}^B_{n+1} = {\bf d}^B_n + \Delta t {\bf v}^B_n,\vspace*{0.25cm}$$ (2.17)

$${\bf v}^A_{n+1} = {\bf v}^A_n +\frac{1}{2} {\Delta t} ({\bf a}^A_n +{\bf a}^A_{n+1}),\vspace*{0.25cm}$$ (2.18)

$${\bf v}^B_{n+1} = {\bf v}^B_n +\frac{1}{2} {\Delta t} ({\bf a}^B_n +{\bf a}^B_{n+1}),\vspace*{0.25cm}$$ (2.19)

$${\bf a}^A_{n+1} = {{\bf M}^A}^{-1}({\bf R}^{A \textrm{in}}_{n+1} + {\bf R}^{A \textrm{co}}_{n+1} + {\bf R}^{A \textrm{ex}}_{n+1} ),\vspace*{0.25cm}$$ (2.20)

$${\bf a}^B_{n+1} = -{\bf a}^B_{n}.\vspace*{0.25cm}$$ (2.21)

At time $t+m \Delta t$, we turn to the explicit update for region B. Using a similar approach presented above, region $B$ is explicitly calculated while region $A$ is approximated through the toggling of its nodal accelerations. According to Smolinski (1989), this update should be performed twice using a time step of $\frac{1}{2} m \Delta t$ or one-half that used in region $A$. This approach was found by Smolinski to correct for the oscillatory nature of the acceleration constraint (toggling) used in the explicit update of region $A$. It is also a reason why the subcycling parameter, $m$, must be an even number. The resulting equations, which are solved twice during the time $t+m \Delta t$, are given by

$${\bf d}^A_{n+1} = {\bf d}^A_n + \frac{1}{2} m \Delta t {\bf v}^A_n,\vspace*{0.25cm}$$ (2.22)

$${\bf d}^B_{n+1} = {\bf d}^B_n + \frac{1}{2} m \Delta t {\bf v}^B_n + \\ \frac{1}{8} m^2 \Delta t^2 {\bf a}^B_n,\vspace*{0.25cm}$$ (2.23)

$${\bf v}^A_{n+1} = {\bf v}^A_n + \frac{1}{4} m {\Delta t} ({\bf a}^A_n + {\bf a}^A_{n+1}),\vspace*{0.25cm}$$ (2.24)

$${\bf v}^B_{n+1} = {\bf v}^B_n + \frac{1}{4} m {\Delta t} ({\bf a}^B_n + {\bf a}^B_{n+1}),\vspace*{0.25cm}$$ (2.25)

$${\bf a}^A_{n+1} = -{\bf a}^A_{n},\vspace*{0.25cm}$$ (2.26)

$${\bf a}^B_{n+1} = {{\bf M}^A}^{-1}({\bf R}^{B \textrm{in}}_{n+1} ... ...bf R}^{B \textrm{co}}_{n+1} + {\bf R}^{B \textrm{ex}}_{n+1} ). \vspace*{0.25cm}$$ (2.27)

As implied above, the bulk of the computational savings associated with the subcycling scheme is achieved when we avoid the calculations of the internal and cohesive force vectors for region $B$ nodes from $t$ to $t+m \Delta t$. These savings may be quite substantial when using more complex constitutive models for the volumetric elements and/or when the number of cohesive elements (and therefore the size of region $A$) is small.

Since the internal force calculations are performed by assembling the local internal force vectors obtained for all the volumetric elements, the most efficient way to take advantage of the subcycling scheme is by flagging the volumetric elements based on the type of nodes ($A$ or $B$) they contain. This allows to quickly ``pass by'' all the elements containing nodes of type $B$ during the first $m$ time steps of the subcycling loop. In the implementation adopted in the present work, all elements of type $A$ have a flag of $1$, all those of type $B$ have a flag of $m$, and all those of ``mixed type'' (i.e., which contain some nodes of type $A$ and some of type $B$) receive a flag of $-1$. This flag then insures that internal force calculations are performed only on those elements having the a flag of $1$ or $-1$, corresponding to elements having at least one node from region $A$. The computational savings we achieve are equal to the number of elements having a flag of $m$, i.e. elements whose internal force calculation is skipped. Figure 2.6 shows an example of applying the element time step flags to a typical mesh.

Figure 2.6: Time step assignment.

An added benefit of the subcycling algorithm is that it can be readily implemented in an adaptive fashion. The time steps of the individual nodes, as well as of the volumetric elements, can be changed at any time during the simulation. This allows the mesh to adapt to the changing critical region. When used in conjunction with adaptive meshing, the time steps are reduced as the local mesh is made finer, and increased as the mesh grows coarser. With dynamic cohesive element insertion, the time steps decrease as cohesive elements are inserted.

Table 2.1: Code implementation of the multi-time step nodal subcycling algorithm (Smolinski, 1989).

Set initial conditions for nodal displacements, velocities and acceleration.

Clear the counter: $c = 0$

Loop over time (conventional explicit central difference method).

if (c $\leq$ m) PRIMARY REGION: $A$ (solve $A$, approximate $B$)

Update $A$ displacements ${\bf d}^A_{n+1}$ using Equation 2.16

Update $B$ displacements ${\bf d}^B_{n+1}$ using Equation 2.17

Calculate ${{\bf R}^A}^\textrm{in}$ and ${\bf R}^\textrm{co}$

Update $A$ accelerations ${\bf a}^A_{n+1}$ using Equation 2.20

Toggle $B$ accelerations ${\bf a}^B_{n+1}$ using Equation 2.21

Update $A$ velocities ${\bf v}^A_{n+1}$ using Equation 2.18

Update $B$ velocities ${\bf v}^B_{n+1}$ using Equation 2.19

Increase counters: $c = c + 1$

Increase time: $t = t + \Delta t$

if (c $>$ m) PRIMARY REGION: $B$ (solve $B$, approximate $A$)

Update $A$ displacements ${\bf d}^A_{n+1}$ using Equation 2.22

Update $B$ displacements ${\bf d}^B_{n+1}$ using Equation 2.23

Calculate ${{\bf R}^B}^\textrm{in}$

Toggle $A$ accelerations ${\bf a}^A_{n+1}$ using Equation 2.26

Update $B$ accelerations ${\bf a}^B_{n+1}$ using Equation 2.27

Update $A$ velocities ${\bf v}^A_{n+1}$ using Equation 2.24

Update $B$ velocities ${\bf v}^B_{n+1}$ using Equation 2.25

Increase counters: $c = c + 1$

If $c = m + 2$, the set $c = 0$