BatheTwoStep

Starting from version 2.7, customisation of spectral radius (${\rho }_{\mathrm{\infty }}$) and sub-step size ($\gamma$) is supported.

The algorithm is known to be able to conserve energy and momentum.

References:

1. 10.1016/j.compstruc.2006.09.004
2. 10.1016/j.compstruc.2018.11.001

The implementation follows the original algorithm. The only difference is the step size defined in suanPan is the mean step size of two sub-step sizes, that is $\mathrm{\Delta }t/2$ in the references.

For further discussions on this type of algorithm, one can check the following references.

1. 10.1007/s11071-023-09065-7

Syntax

integrator BatheTwoStep (1) [2] [3]
# (1) int, unique integrator tag
# [2] double, spectral radius, \rho_\infty, default: 0
# [2] double, sub-step size, gamma, default: 0.5

Using integrator BatheTwoStep (1) with two optional parameters omitted gives the same results as in versions prior to version 2.7.

The spectral radius ${\rho }_{\mathrm{\infty }}$ ranges from 0 to 1, both ends inclusive.

The sub-step size $\gamma$ ranges from 0 to 1, both ends exclusive.

Theory

The First Sub-step

For trapezoidal rule,

$$v_{n+1}=v_n+{\displaystyle \frac{\gamma \mathrm{\Delta }t}{2}}(a_n+a_{n+1}),$$

$$u_{n+1}=u_n+{\displaystyle \frac{\gamma \mathrm{\Delta }t}{2}}(v_n+v_{n+1}).$$

Then,

$$u_{n+1}=u_n+{\displaystyle \frac{\gamma \mathrm{\Delta }t}{2}}(v_n+v_n+{\displaystyle \frac{\gamma \mathrm{\Delta }t}{2}}(a_n+a_{n+1})),$$

$$\mathrm{\Delta }u=\gamma \mathrm{\Delta }t{v}_n+{\displaystyle \frac{{\gamma }^2\mathrm{\Delta }{t}^2}{4}}{a}_n+{\displaystyle \frac{{\gamma }^2\mathrm{\Delta }{t}^2}{4}}{a}_{n+1}.$$

One could obtain

$$a_{n+1}={\displaystyle \frac{4}{{\gamma }^2\mathrm{\Delta }{t}^2}}\mathrm{\Delta }u-{\displaystyle \frac{4}{\gamma \mathrm{\Delta }t}}{v}_n-a_n,\mathrm{\Delta }a={\displaystyle \frac{4}{{\gamma }^2\mathrm{\Delta }{t}^2}}\mathrm{\Delta }u-{\displaystyle \frac{4}{\gamma \mathrm{\Delta }t}}{v}_n-2a_n,$$

$$v_{n+1}={\displaystyle \frac{2}{\gamma \mathrm{\Delta }t}}\mathrm{\Delta }u-v_n,\mathrm{\Delta }v={\displaystyle \frac{2}{\gamma \mathrm{\Delta }t}}\mathrm{\Delta }u-2v_n.$$

The effective stiffness is then

$$\overline{K}=K+{\displaystyle \frac{2}{\gamma \mathrm{\Delta }t}}C+{\displaystyle \frac{4}{{\gamma }^2\mathrm{\Delta }{t}^2}}M.$$

The Second Sub-step

The second step is computed by

$$v_{n+2}=v_n+\mathrm{\Delta }t(q_0{a}_n+q_1{a}_{n+1}+q_2{a}_{n+2}),$$

$$u_{n+2}=u_n+\mathrm{\Delta }t(q_0{v}_n+q_1{v}_{n+1}+q_2{v}_{n+2}).$$

The parameters satisfy $q_0+q_1+q_2=1$, and

$$q_1={\displaystyle \frac{{\rho }_{\mathrm{\infty }}+1}{2\gamma ({\rho }_{\mathrm{\infty }}-1)+4}},q_2=0.5-\gamma q_1.$$

Hence,

$$a_{n+2}={\displaystyle \frac{v_{n+2}-v_n}{q_2\mathrm{\Delta }t}}-{\displaystyle \frac{q_0}{q_2}}{a}_n-{\displaystyle \frac{q_1}{q_2}}{a}_{n+1},$$

$$v_{n+2}={\displaystyle \frac{u_{n+2}-u_n}{q_2\mathrm{\Delta }t}}-{\displaystyle \frac{q_0}{q_2}}{v}_n-{\displaystyle \frac{q_1}{q_2}}{v}_{n+1}.$$

The effective stiffness is then

$$\overline{K}=K+{\displaystyle \frac{1}{q_2\mathrm{\Delta }t}}C+{\displaystyle \frac{1}{q_2^2\mathrm{\Delta }{t}^2}}M.$$