Composite methods for structural dynamics

From HandWiki

Composite methods are an approach applied in structural dynamics and related fields. They combine various methods in each time step, in order to acquire the advantages of different methods. The existing composite methods show satisfactory accuracy and powerful numerical dissipation, which is particularly useful for solving stiff problems[1] and differential-algebraic equations.[2]

Definitions

After spatial discretization, structural dynamics problems are generally described by the second-order ordinary differential equation:

[math]\displaystyle{ M\ddot{u} + C\dot{u} + f(u,t) = R(t) }[/math].

Here [math]\displaystyle{ u }[/math], [math]\displaystyle{ \dot{u} }[/math] and [math]\displaystyle{ \ddot{u} }[/math] denote the displacement, velocity and acceleration vectors respectively, [math]\displaystyle{ M }[/math] is the mass matrix, [math]\displaystyle{ C }[/math] is the damping matrix, [math]\displaystyle{ f(u,t) }[/math] collects the internal force, and [math]\displaystyle{ R(t) }[/math] is the external load. At the initial time [math]\displaystyle{ t_0 }[/math], the initial displacement and velocity are supposed to be given as [math]\displaystyle{ u_0 }[/math] and [math]\displaystyle{ \dot{u}_0 }[/math], respectively, and the initial acceleration can be solved as

[math]\displaystyle{ \ddot{u}_0=M^{-1}(R(t_0)-C\dot{u}_0-f(u_0,t_0)) }[/math].

For numerical analysis, the overall time domain [math]\displaystyle{ [t_0,t_N] }[/math] is divided into a series of time steps by [math]\displaystyle{ t_1 }[/math], [math]\displaystyle{ t_2 }[/math], [math]\displaystyle{ \cdots }[/math], [math]\displaystyle{ t_k }[/math], [math]\displaystyle{ t_{k+1} }[/math], [math]\displaystyle{ \cdots }[/math]. Taking the step [math]\displaystyle{ [t_k,t_{k+1}] }[/math] ([math]\displaystyle{ t_{k+1}-t_k=h }[/math] is the step size), the main concept of composite methods is to subdivide the current step to several sub-steps [math]\displaystyle{ [t_k,t_{k+\gamma_1}] }[/math], [math]\displaystyle{ [t_{k+\gamma_1},t_{k+\gamma_2}] }[/math], [math]\displaystyle{ \cdots }[/math], and to use different numerical methods in each sub-step.

Although there are lots of available methods, see the review,[3] the existing composite methods basically employ the combination of the trapezoidal rule and linear multistep methods. However, to acquire at least second-order accuracy and unconditional stability, the scalar parameters of each method and the division of sub-steps need to be determined carefully.

Two examples of composite method

Two-sub-step Bathe method

The Bathe method [4][5] is a two-sub-step method. In the first sub-step [math]\displaystyle{ [t_k,t_{k+\gamma}] }[/math] ([math]\displaystyle{ t_{k+\gamma}-t_k=\gamma h }[/math], [math]\displaystyle{ \gamma\in(0,1) }[/math]), the trapezoidal rule is used as:

[math]\displaystyle{ u_{k+\gamma}=u_{k}+\frac{\gamma h}{2}(\dot{u}_k+\dot{u}_{k+\gamma}) }[/math]

[math]\displaystyle{ \dot{u}_{k+\gamma}=\dot{u}_{k}+\frac{\gamma h}{2}(\ddot{u}_k+\ddot{u}_{k+\gamma}) }[/math]

[math]\displaystyle{ M\ddot{u}_{k+\gamma}+C\dot{u}_{k+\gamma}+f(u_{k+\gamma},t_{k+\gamma})=R(t_{k+\gamma}) }[/math]

In the second sub-step [math]\displaystyle{ [t_{k+\gamma},t_{k+1}] }[/math] ([math]\displaystyle{ t_{k+1}-t_{k+\gamma}=(1-\gamma) h }[/math]), the 3-point Euler backward method is employed as

[math]\displaystyle{ \dot{u}_{k+1}=\frac{1-\gamma}{\gamma h}u_k-\frac{1}{(1-\gamma)\gamma h}u_{k+\gamma}+\frac{2-\gamma}{(1-\gamma)h}u_{k+1} }[/math]

[math]\displaystyle{ \ddot{u}_{k+1}=\frac{1-\gamma}{\gamma h}\dot{u}_k-\frac{1}{(1-\gamma)\gamma h}\dot{u}_{k+\gamma}+\frac{2-\gamma}{(1-\gamma)h}\dot{u}_{k+1} }[/math]

[math]\displaystyle{ M\ddot{u}_{k+1}+C\dot{u}_{k+1}+f(u_{k+1},t_{k+1})=R(t_{k+1}) }[/math]

For nonlinear dynamics, that is, the internal force [math]\displaystyle{ f }[/math] is a nonlinear function with respect to [math]\displaystyle{ u }[/math], the Newton-Raphson iterations can be used to solve the nonlinear equations per step. The parameter [math]\displaystyle{ \gamma }[/math] is usually set as [math]\displaystyle{ \frac{1}{2} }[/math] and [math]\displaystyle{ 2-\sqrt{2} }[/math] in practice.

The Bathe method is second-order accurate and unconditionally stable from linear analysis. Besides, this method can provide strong numerical dissipation for high-frequency content, which is helpful to damp out the stiff components and enhance the stability for nonlinear dynamics.

On this basis, to acquire prescribed degree of numerical dissipation, the [math]\displaystyle{ \rho_\infty }[/math]-Bathe method [6] was developed by replacing the 3-point Euler backward method in the second sub-step with a general formula:

[math]\displaystyle{ u_{k+1}=u_k+h(q_0\dot{u}_k+q_1\dot{u}_{k+\gamma}+q_2\dot{u}_{k+1}) }[/math]

[math]\displaystyle{ \dot{u}_{k+1}=\dot{u}_k+h(q_0\ddot{u}_k+q_1\ddot{u}_{k+\gamma}+q_2\ddot{u}_{k+1}) }[/math]

[math]\displaystyle{ M\ddot{u}_{k+1}+C\dot{u}_{k+1}+f(u_{k+1},t_{k+1})=R(t_{k+1}) }[/math]

The parameters are selected as recommended

[math]\displaystyle{ \gamma=\frac{2-\sqrt{2(1+\rho_\infty)}}{1-\rho_\infty}\text{ if }\rho_\infty\in[0,1); \gamma=\frac{1}{2} \text{ if } \rho_\infty=1 }[/math]

[math]\displaystyle{ q_1=\frac{\rho_\infty+1}{2\gamma(\rho_\infty-1)+4},q_0=(\gamma-1)q_1+\frac{1}{2},q_2=-\gamma q_1+\frac{1}{2} }[/math]

With the set of parameters, the [math]\displaystyle{ \rho_\infty }[/math]-Bathe method can also achieve second-order accuracy and unconditional stability. Moreover, by adjusting the parameter [math]\displaystyle{ \rho_\infty }[/math], this method can provide tunable degree of numerical dissipation. The method with a smaller [math]\displaystyle{ \rho_\infty }[/math] shows stronger numerical dissipation, but lower accuracy in the low-frequency content. When [math]\displaystyle{ \rho_\infty=0 }[/math], it is equivalent to the original Bathe method with [math]\displaystyle{ \gamma=2-\sqrt{2} }[/math].

Three-sub-step composite method

Following the idea of the Bathe method, the three-sub-step composite methods that use the trapezoidal rule in the first two sub-steps were also discussed.[7][8][9] They divides the current step into [math]\displaystyle{ [t_k,t_{k+\gamma_1}] }[/math], [math]\displaystyle{ [t_{k+\gamma_1},t_{k+\gamma_2}] }[/math] and [math]\displaystyle{ [t_{k+\gamma_2},t_{k+1}] }[/math], and generally, the first two sub-steps are set as equal size, that is [math]\displaystyle{ \gamma_2=2\gamma_1 }[/math]. In the first two sub-steps, the trapezoidal rule is used, as

[math]\displaystyle{ u_{k+\gamma_1}=u_{k}+\frac{\gamma_1 h}{2}(\dot{u}_k+\dot{u}_{k+\gamma_1}) }[/math]

[math]\displaystyle{ \dot{u}_{k+\gamma_1}=\dot{u}_{k}+\frac{\gamma_1 h}{2}(\ddot{u}_k+\ddot{u}_{k+\gamma_1}) }[/math]

[math]\displaystyle{ M\ddot{u}_{k+\gamma_1}+C\dot{u}_{k+\gamma_1}+f(u_{k+\gamma_1},t_{k+\gamma_1})=R(t_{k+\gamma_1}) }[/math]

and

[math]\displaystyle{ u_{k+\gamma_2}=u_{k+\gamma_1}+\frac{(\gamma_2-\gamma_1) h}{2}(\dot{u}_{k+\gamma_1}+\dot{u}_{k+\gamma_2}) }[/math]

[math]\displaystyle{ \dot{u}_{k+\gamma_2}=\dot{u}_{k+\gamma_1}+\frac{(\gamma_2-\gamma_1) h}{2}(\ddot{u}_{k+\gamma_1}+\ddot{u}_{k+\gamma_2}) }[/math]

[math]\displaystyle{ M\ddot{u}_{k+\gamma_2}+C\dot{u}_{k+\gamma_2}+f(u_{k+\gamma_2},t_{k+\gamma_2})=R(t_{k+\gamma_2}) }[/math]

In the last sub-step, a general formula is utilized as

[math]\displaystyle{ u_{k+1}=u_k+h(c_0\dot{u}_k+c_1\dot{u}_{k+\gamma_1}+c_2\dot{u}_{k+\gamma_2}+c_3\dot{u}_{k+1}) }[/math]

[math]\displaystyle{ \dot{u}_{k+1}=\dot{u}_k+h(c_0\ddot{u}_k+c_1\ddot{u}_{k+\gamma_1}+c_2\ddot{u}_{k+\gamma_2}+c_3\ddot{u}_{k+1}) }[/math]

[math]\displaystyle{ M\ddot{u}_{k+1}+C\dot{u}_{k+1}+f(u_{k+1},t_{k+1})=R(t_{k+1}) }[/math]

For this method, Li et al.[8] offered two optimal set of parameters, as

[math]\displaystyle{ a=\frac{1}{2}(1\pm\rho_\infty),c_0=\frac{-(a+1)\gamma_1^2+4\gamma_1-1}{4\gamma_1},c_1=\frac{1+(a-1)\gamma_1}{2},c_2=\frac{(1-a)\gamma_1^2-2\gamma_1+1}{4\gamma_1},c_3=\frac{\gamma_1}{2} }[/math]

Here [math]\displaystyle{ \gamma_2=2\gamma_1 }[/math] is assumed, and [math]\displaystyle{ \gamma_1 }[/math] is the minimum value that satisfies [math]\displaystyle{ 48(a-1)\gamma_1^4-32(a-5)\gamma_1^3-192\gamma_1^2+96\gamma_1-16\geq0 }[/math].

The resulting two sub-families are all second-order accurate, unconditionally stable, and can provide tunable numerical dissipation by adjusting [math]\displaystyle{ \rho_\infty }[/math]. They become the same when [math]\displaystyle{ \rho_\infty=0 }[/math]. When [math]\displaystyle{ 0\lt \rho_\infty\lt 1 }[/math], the sub-family with [math]\displaystyle{ a=\frac{1}{2}(1-\rho_\infty) }[/math] shows better amplitude and period accuracy than the [math]\displaystyle{ \rho_\infty }[/math]-Bathe method under the same computational costs, and the sub-family with [math]\displaystyle{ a=\frac{1}{2}(1+\rho_\infty) }[/math] further improves the period accuracy at the cost of lower amplitude accuracy.

Analysis

In structural dynamics, the test model for property analysis is the single degree-of-freedom homogeneous equation, as

[math]\displaystyle{ \ddot{u}+2\xi\omega\dot{u}+\omega^2u=0 }[/math]

Here [math]\displaystyle{ \xi }[/math] is the damping ratio and [math]\displaystyle{ \omega }[/math] is the natural frequency. Applying the composite method to the test model yields the compact scheme

[math]\displaystyle{ X_{k+1}=AX_{k} }[/math]

Here [math]\displaystyle{ X_{k}=\{u_k;\dot{u}_k;\ddot{u}_k\} }[/math] and [math]\displaystyle{ A }[/math] is the amplitude matrix, which governs the properties of a method. Generally, [math]\displaystyle{ A }[/math] has one zero characteristic root and a pair of conjugate complex roots [math]\displaystyle{ \lambda_{1,2} }[/math], which can be solved from

[math]\displaystyle{ \lambda^2-A_1\lambda+A_2=0 }[/math]

Here [math]\displaystyle{ A_1 }[/math] is the trace of [math]\displaystyle{ A }[/math] and [math]\displaystyle{ A_2 }[/math] is the sum of second-order principal minors of [math]\displaystyle{ A }[/math]. They are functions of [math]\displaystyle{ \xi }[/math], [math]\displaystyle{ \omega h }[/math], and the parameters of the method.

Accuracy

From the compact scheme, the difference equation only with respect to the displacement can be written as

[math]\displaystyle{ u_{k+1}-A_1u_{k}+A_2u_{k-1}=0 }[/math]

The local truncation error [math]\displaystyle{ \sigma }[/math] is defined as

[math]\displaystyle{ \sigma=u(t_{k+1})-A_1u(t_{k})+A_2u(t_{k-1}) }[/math]

The method is called [math]\displaystyle{ s }[/math]th-order accurate if [math]\displaystyle{ \sigma=O(h^{s+1}) }[/math].

Stability

For physically stable systems ([math]\displaystyle{ \xi\geq0 }[/math], [math]\displaystyle{ \omega\geq0 }[/math]), the method can give stable solutions if the spectral radius [math]\displaystyle{ \rho=\max\{|\lambda|\}\leq1 }[/math]. A method is called unconditionally stable if the condition [math]\displaystyle{ \rho\leq1 }[/math] is satisfied for any [math]\displaystyle{ h\geq0 }[/math], otherwise it is called conditionally stable. The spectral radius at the high-frequency limit, i.e. [math]\displaystyle{ \omega h\rightarrow+\infty }[/math], is denoted as [math]\displaystyle{ \rho_\infty }[/math], which is usually employed to indicate the degree of numerical dissipation, as used above.

Amplitude decay ratio and period elongation ratio

In addition to the accuracy order, the amplitude decay ratio and period elongation ratio are also usually evaluated to measure the amplitude and period accuracy in the low-frequency content. The exact solution of the test model is

[math]\displaystyle{ u(t)=\text{e}^{-\xi\omega t}(c_1\cos\omega_dt+c_2\sin\omega_dt),\omega_d=\omega\sqrt{1-\xi^2} }[/math]

Here [math]\displaystyle{ c_1 }[/math] and [math]\displaystyle{ c_2 }[/math] are constants determined by the initial conditions. The numerical solution can be also expressed as a similar form, as

[math]\displaystyle{ u_k=\text{e}^{-\overline{\xi}\overline{\omega}t_k}(\overline{c}_1\cos\overline{\omega}_dt_k+\overline{c}_2\sin\overline{\omega}_dt_k),\overline{\omega}_d=\overline{\omega}\sqrt{1-\overline{\xi}^2} }[/math]

Likewise, [math]\displaystyle{ \overline{c}_1 }[/math] and [math]\displaystyle{ \overline{c}_2 }[/math] are also determined by the initial conditions and they should be close to [math]\displaystyle{ c_1 }[/math] and [math]\displaystyle{ c_2 }[/math] respectively for a convergent method. The damping ratio [math]\displaystyle{ \overline{\xi} }[/math] and frequency [math]\displaystyle{ \overline{\omega} }[/math] can be obtained from the norm [math]\displaystyle{ |\lambda| }[/math] and phase [math]\displaystyle{ \angle\lambda }[/math], as[10]

[math]\displaystyle{ \overline{\xi}=-\frac{\ln|\lambda|}{\sqrt{(\angle\lambda)^2+(\ln|\lambda|)^2}},\overline{\omega}=\frac{\sqrt{(\angle\lambda)^2+(\ln|\lambda|)^2}}{h} }[/math]

Here [math]\displaystyle{ \overline{\xi} }[/math] is called the amplitude decay ratio, and [math]\displaystyle{ \frac{\overline{T}-T}{T} }[/math] ([math]\displaystyle{ \overline{T}=\frac{2\pi}{\overline{\omega}},T=\frac{2\pi}{\omega} }[/math]) is called the period elongation ratio.

Example

File:Spectral radius of the Bathe method.tif

File:Amplitude decay ratio of the Bathe method.tif

File:Period elongation ratio of the Bathe method.tif

Consider the Bathe method, [math]\displaystyle{ A_1 }[/math] and [math]\displaystyle{ A_2 }[/math] have the form as

[math]\displaystyle{ A_1=\frac{2(\gamma^4-4\gamma^3+6\gamma^2-4)\omega^2h^2+8(\gamma-2)^2}{(\gamma^2\omega^2h^2+4)((\gamma-1)^2\omega^2h^2+(\gamma-2)^2)} }[/math]

[math]\displaystyle{ A_2=\frac{(\gamma^4-4\gamma^3+8\gamma^2-8\gamma+4)\omega^2h^2+4(\gamma-2)^2}{(\gamma^2\omega^2h^2+4)((\gamma-1)^2\omega^2h^2+(\gamma-2)^2)} }[/math]

Here the undamped case, i.e. [math]\displaystyle{ \xi=0 }[/math], is considered for simplicity. One can check that this method can satisfy the conditions of second-order accuracy and unconditional stability. With [math]\displaystyle{ \gamma=\frac{1}{2} }[/math] and [math]\displaystyle{ 2-\sqrt{2} }[/math], the spectral radius, amplitude decay ratio, and period elongation ratio are shown here. It can be observed that this method can provide good amplitude and period accuracy in the low-frequency content, while strong numerical dissipation, as [math]\displaystyle{ \rho_\infty=0 }[/math], in the high-frequency content.

See also

References

  1. ↑ Hairer, Ernst; Wanner, Gerhard (1996). Solving ordinary differential equations (Second ed.). Berlin: Springer-Verlag. ISBN 978-3-540-60452-5. 
  2. ↑ Kunkel, Peter; Mehrmann, Volker Ludwig (2006). Differential-algebraic equations : analysis and numerical solution. European Mathematical Society. ISBN 978-3-03719-017-3. 
  3. ↑ Tamma, Kumar K.; Har, Jason; Zhou, Xiangmin; Shimada, Masao; Hoitink, Andrew (15 July 2011). "An Overview and Recent Advances in Vector and Scalar Formalisms: Space/Time Discretizations in Computational Dynamics—A Unified Approach". Archives of Computational Methods in Engineering 18 (2): 119–283. doi:10.1007/s11831-011-9060-y. 
  4. ↑ Bathe, Klaus-Jürgen; Baig, Mirza M. Irfan (December 2005). "On a composite implicit time integration procedure for nonlinear dynamics". Computers & Structures 83 (31–32): 2513–2524. doi:10.1016/j.compstruc.2005.08.001. 
  5. ↑ Bathe, Klaus-Jürgen; Noh, Gunwoo (May 2012). "Insight into an implicit time integration scheme for structural dynamics". Computers & Structures 98-99: 1–6. doi:10.1016/j.compstruc.2012.01.009. 
  6. ↑ Noh, Gunwoo; Bathe, Klaus-Jürgen (February 2019). "The Bathe time integration method with controllable spectral radius: The ρ∞-Bathe method". Computers & Structures 212: 299–310. doi:10.1016/j.compstruc.2018.11.001. 
  7. ↑ Chandra, Yenny; Zhou, Yang; Stanciulescu, Ilinca; Eason, Thomas; Spottswood, Stephen (11 April 2015). "A robust composite time integration scheme for snap-through problems". Computational Mechanics 55 (5): 1041–1056. doi:10.1007/s00466-015-1152-3. Bibcode: 2015CompM..55.1041C. 
  8. ↑ 8.0 8.1 Li, Jinze; Yu, Kaiping; Li, Xiangyang (26 April 2019). "A novel family of controllably dissipative composite integration algorithms for structural dynamic analysis". Nonlinear Dynamics 96 (4): 2475–2507. doi:10.1007/s11071-019-04936-4. 
  9. ↑ Ji, Yi; Xing, Yufeng (April 2020). "An optimized three-sub-step composite time integration method with controllable numerical dissipation". Computers & Structures 231: 106210. doi:10.1016/j.compstruc.2020.106210. 
  10. ↑ Zhou, X.; Tamma, K. K. (7 February 2004). "Design, analysis, and synthesis of generalized single step single solve and optimal algorithms for structural dynamics". International Journal for Numerical Methods in Engineering 59 (5): 597–668. doi:10.1002/nme.873. Bibcode: 2004IJNME..59..597Z. 



Retrieved from "https://handwiki.org/wiki/index.php?title=Composite_methods_for_structural_dynamics&oldid=2628261"