Numerical Integration Methods for Multibody Dynamics

1 Numerical Integration Methods for Multibody Dynamics Lesson 8. Numerical Integration Methods for Multibody Dynamics Accuracy of Explicit Euler Method By Taylor series expansions, Explicit Euler Method 2 3. +1 = + + + . 2 6 . 2. +1 = + + + Truncation error 2. +1 = + + Truncation error 3. Truncation error ( ) < ( : user defined error tolerance). 6. 3 6 . Allowable step size .. must be estimated). ( . 3. Accuracy of Implicit Euler Method By Taylor series expansions, Implicit Euler Method 2 3. +1 = + + + . 2 6 . 2 2 3. +1 = + + +1 + Truncation error = + + +1 + . 2 2 3. where, +1 = + +1 + Truncation error 2 . = =. 2 . + +. 2. + . 2 +1 +1. 2 . 2 . 2 . = +1 + . 2. 3. Truncation error ( ) < ( : user defined error tolerance). 3. 3 3 . Allowable step size .. must be estimated). ( . 4. Stability of Euler Methods . Explicit Method.

2 +1 = + = . 1. +1 = (1 + ) |1 + | < 1. 1 = 1 + 0 = (1 + ) 0 inside stable Implicit Method +1 = + +1 = . 1 1. +1 = | | < 1 (Absolutely stable). 1 1 . 1 1. 1 = = . 1 0 (1 ) 0. 5. Numerical Damping of a Numerical Integrator Why Numerical damping ? To provide Numerical stability without sacrificing accuracy for discontinuous or chaotic systems Numerical damping of Backward Euler method +1 = + +1 = . 1 1. +1 = | | < 1 (Absolutely stable). 1 1 . 1 1. 1 = = . 1 0 (1 ) 0. Bigger step size provides more Numerical stability Therefore, It is important to have a small step size not to disturb accuracy with good Numerical stability 6. Effect of Step Size and Numerical Damping MK System m = 1 kg, k = 100 N/mm 0 = 100 mm, 0 = 10000 mm/s Ansys Motion Results Numerical Damping = Numerical Damping = Numerical Damping = Position (mm).

3 Position (mm). Position (mm). Time (s) Time (s) Time (s). Bigger step size & greater Numerical damping provides more Numerical stability Therefore, It is important to have a small step size not to disturb accuracy with good Numerical stability 7. Numerical Integration Algorithm of an Example . 3 3 Ground and body A are connected by a translational joint 1 3 Body A and body B are connected by a revolute joint C Body B and body C are connected by a beam element 1 . 2 2 . 0 1 2 3.. 2. 1 = 0 : 1 = 1 : 2 = 2 : 3 = 3. 1 . 0 1 2 3.. 0. Constraint equation: rev = + + = . 1 cos 1 sin 1 0 2 cos 2 sin 2 0. = + 2 +. 1 sin 1 cos 1 0 sin 2 cos 2 0. 1 2. = = . 1 2. 8. Constraints and Their Jacobian Constraint equation (continued). 1 0 cos 1 sin 1 0. T. ( 1 ( 1 ) 0 ) 1 ( 1 ) 1 1 0 sin 1. 0 1 sin 1 cos 1 1. trans = T = = = . 01 1 ( 1 ) 1 1 1 cos 1 sin 1 0 1 cos 1 1 sin 1.

4 Sin 1 cos 1 1. 1 2. 1 2. = sin 1 = . 1 cos 1 1 sin 1. Jacobian Matrix T T. General coordinates : = 1 T 2 T 3 T = 1 1 1 2 2 2 3 3 3. 1 0 0 1 0 0 0 0 0. 0 1 0 0 1 0 0 0 0. =. 0 0 cos 1 0 0 0 0 0 0. sin 1 cos 1 x1 cos 1 y1 sin 1 0 0 0 0 0 0. 9. 9. Generalized Forces Beam force 2 = 2,t+ t 2,t 23 = 2 T 3 2 02. 3 = 3,t+ t 3,t Y. 3 Beam = 2 23 = ( 2 2 T 3 2 02 ). 3. Beam = [( 2 2 T 3 2 02 )]. 2 [ G 3 2 02 ] ( 2 2 T = G ). 2 2 . 2 = G 2 G 3. 2 02. 2 2. = [ Beam ] = G G . 3 3. X. 12 . 11 0 0 11 = , 22 =. 3. , = 0 22 26. 6 4 . 0 26 66 26 = 62 = 2. , 66 =.. Spring force 01 1 0 T. 1 0. spring = = 1 0 =. 01 x1 0. 10. Equations of Motion Equations of Motion T . =.. T . 1 0 0 0 0 0 0 0 0 1 0 0 sin 1 1 ( 1 0 ). 0 1 0 0 0 0 0 0 0 0 1 0 cos 1 1 1 . 0 0 1 0 0 0 0 0 0 0 0 cos 1 1 sin 1 1 cos 1 1 0. 0 0 0 2 0 0 0 0 0 1 0 0 0 2 , . 0 0 0 0 2 0 0 0 0 0 1 0 0 2 2 +.

5 0 0 0 0 0 2 0 0 0 0 0 0 0 2 , . 0 0 0 0 0 0 3 0 0 0 0 0 0 3 = , . 0 0 0 0 0 0 0 3 0 0 0 0 0 3 3 , . 0 0 0 0 0 0 0 0 3 0 0 0 0 3 , . 1 0 0 1 0 0 0 0 0 0 0 0 0 1 1. 0 1 0 0 1 0 0 0 0 0 0 0 0 2 2. 0 0 cos 1 0 0 0 0 0 0 0 0 0 0 3 3. sin 1 cos 1 1 cos 1 1 sin 1 0 0 0 0 0 0 0 0 0 0 4 4.. = = . T . = = + = .. = . 11. Numerical Integration Algorithm by Explicit Method Assume step time size : h = At t=0 , set initial values 0 and 0.. Solve = and find 0 and using 0 and 0.. Numerical Integration 2. 1 = 0 + 0 + 0. 2. 1 = 0 + 0. Increase time by Repeat - . 12. Initialization of Implicit Method Initialization of Implicit method Solve 0 from EQM with given initial 0 and 0.. =.. Update = + . Set initial guesses +1, +1, +1 and +1. Calculate residuals ~ . 1 = +1 +1. 2 = +1 +1. 3 = +1 + T +1 +1. 4 = +1. 13. Jacobian Matrix for Newton's Method Find new +1 , +1 , +1 and +1 by Newton's Method 1 1 1 1.

6 2 2 2 2.. = . = 3 3 3 3. = T. T.. 4 4 4 4.. T T. where = 1 2 3 4 , = . =. 14. Newton's Method Iteration Perform Newton's method iteration Solve the matrix equation ( ) = . +1 ( ) ( ) +1 ( ). +1 ( ) ( ) +1 ( ). T. T = . +1 ( ) + T +1 ( ) +1 ( ). +1 ( ). Update +1 ( +1) , +1 ( +1), +1 ( +1) and +1. +1 . +1 . = +. +1 . +1 . Repeat this step until converge to the exact solution within a given tolerance Repeat - . 15. Solver Options for Ansys Motion Auto Default option, the solver type is automatically changed by a given condition. Frontal solver The Nested dissection method is used to eliminate intermediate variables. This solver has shown a better performance in HPC machine than Sparse solver. This solver generally uses a smaller memory than Sparse solver. Super solver This solver uses the frontal solver. LU decomposition values are saved and reused as long as Newton iteration is converged.

7 This solver performs well for slowly moving contact or small deformation problems. In current HPC implementation, each flexible body is assigned to different core. Therefore, it performs well for a problem with many flexible bodie. 16. Solver Options for Ansys Motion Sparse solver A matrix is represented by the Compressed Sparse Row(CSR) or Compressed Row Storage(CRS) format. Intel MKL library, a commercial library, is used. SMP version performs similar as frontal solver. PCG. Preconditioned Conjugate Gradient (PCG) iterative solver. This solver performs well for small deformation models. This solver generally takes longer than general linear solver for shell element or large deformation problems. Preprocessing time is required to find a good preconditioner. 17.