Research Papers

A Matrix-Free Newton–Krylov Parallel Implicit Implementation of the Absolute Nodal Coordinate Formulation

[+] Author and Article Information
Daniel Melanz, Naresh Khude, Dan Negrut

Department of Mechanical Engineering,
1513 University Avenue,
Madison, WI 53706

Paramsothy Jayakumar

U.S. Army Tank Automotive Research,
Development and Engineering Center,
6501 East 11 Mile Road,
Warren, MI 48397

Contributed by the Design Engineering Division of ASME for publication in the JOURNAL OF COMPUTATIONAL AND NONLINEAR DYNAMICS. Manuscript received December 3, 2012; final manuscript received May 1, 2013; published online October 9, 2013. Assoc. Editor: Hiroyuki Sugiyama.

J. Comput. Nonlinear Dynam 9(1), 011006 (Oct 09, 2013) (9 pages) Paper No: CND-12-1214; doi: 10.1115/1.4025281 History: Received December 03, 2012; Revised May 01, 2013

This paper sets out to demonstrate three things: (i) implicit integration with absolute nodal coordinate formulation (ANCF) is effective in handling very stiff systems when an accurate computation of the sensitivity matrix is part of the solution sequence, (ii) parallel computing can provide a vehicle for ANCF to tackle very large kinematically constrained problems with millions of degrees of freedom and produce results in a matter of seconds, and (iii) large systems of equations associated with implicit integration can be solved in parallel by relying on an iterative approach that avoids costly matrix factorizations, which would be prohibitively expensive and memory intensive. For (iii), the approach adopted relies on a Krylov–subspace method that is invoked in the Newton stage at each time step of the numerical solution process. The proposed approach is validated against a commercial package and several simple systems for which analytical solutions are available. A set of numerical experiments demonstrates the scaling of the parallel solution method and provides insights in relation to the size of ANCF problems that are tractable using graphics processing unit (GPU) parallel computing and implicit numerical integration.

Copyright © 2014 by ASME
Your Session has timed out. Please sign back in to continue.


Von Dombrowski, S., 2002, “Analysis of Large Flexible Body Deformation in Multibody Systems Using Absolute Coordinates,” Multibody Syst. Dyn., 8, pp. 409–432. [CrossRef]
Dufva, K., and Shabana, A. A., 2005, “Analysis of Thin Plate Structures Using the Absolute Nodal Coordinate Formulation,” Proc. IMechE Part K: J Multibody Dyn., 219, pp. 345–355.
Gerstmayr, J., and Shabana, A., 2006, “Analysis of Thin Beams and Cables Using the Absolute Nodal Co-Ordinate Formulation,” Nonlinear Dyn., 45, pp. 109–130. [CrossRef]
Shabana, A. A., 2008, Computational Continuum Mechanics, Cambridge University Press, New York.
Melanz, D., 2012, “On the Validation and Applications of a Parallel Flexible Multi-body Dynamics Implementation,” M.S. thesis, University of Wisconsin-Madison, Madison, WI.
Shabana, A. A., 2005, Dynamics of Multibody Systems, 3rd ed., Cambridge University Press, New York.
Haug, E. J., 1989, Computer-Aided Kinematics and Dynamics of Mechanical Systems. Volume I: Basic Methods, Allyn and Bacon, Boston, MA.
Haug, E. J., Wu, S., and Yang, S., 1986, “Dynamics of Mechanical Systems With Coulomb Friction, Stiction, Impact, and Constraint Addition-Deletion-I,” Mech. Mach. Theory, 21, pp. 401–406. [CrossRef]
Hairer, E., and Wanner, G., 1996, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Vol. 14, 2nd ed., Springer, Berlin.
Khude, N., Heyn, T., Jay, L.O., and Negrut, D., 2007, “A Comparison of Low Order Numerical Integration Formulas for Time Domain Analysis of Mechanical Systems,” Proceedings of the ECCOMAS Thematic Conference, Milan, Italy.
Negrut., D., Rampalli, R., Ottarsson, G., and Sajdak, A., 2007, “On an Implementation of the HHT Method in the Context of Index 3 Differential Algebraic Equations of Multibody Dynamics,” ASME J. Comput. Nonlinear Dyn., 2(1), pp. 73–85. [CrossRef]
Hussein, B., Negrut, D., and Shabana, A., 2008, “Implicit and Explicit Integration in the Solution of the Absolute Nodal Coordinate Differential/Algebraic Equations,” Nonlinear Dyn., 54, pp. 283–296. [CrossRef]
Newmark, N. M., 1959, “A Method of Computation for Structural Dynamics,” ASCE J. Eng. Mech., pp. 67–94.
Arnold, M., and Bruls, O., 2007, “Convergence of the Generalized-Alpha Scheme for Constrained Mechanical Systems,” Multibody Syst. Dyn., 18, pp. 185–202. [CrossRef]
Khude, N., Jay, L., Schaffer, A., and Negrut, D., 2007, “A Discussion of Low Order Numerical integration Formulas for Rigid and Flexible Multibody Dynamics,” Proceeding of the 6th ASME International Conference on Multibody Systems, Nonlinear Dynamics and Control, Las Vegas, NV, Paper No. DETC2007-35666.
Bottasso, C. L., Bauchau, O., and Cardona, A., 2007, “Time-Step-Size-Independent Conditioning and Sensitivity to Perturbations in the Numerical Solution of Index Three Differential Algebraic Equations,” Siam J. Sci. Comput., 29, pp. 397–414. [CrossRef]
van derVorst, H. A., 1992, “BI-CGSTAB: a Fast and Smoothly Converging Variant of BI-CG for the Solution of Nonsymmetric Linear Systems,” SIAM J. Sci. Stat. Comput., 13, pp. 631–644. [CrossRef]
Atkinson, K. E., 1989, An Introduction to Numerical Analysis, 2nd ed., Wiley, New York.
Saad, Y., 2003, Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics, Philadelphia, PA.
NVIDIA, 2011, “Compute Unified Device Architecture Programming Guide 4.0,” Available at http://developer.download.nvidia.com/compute/cuda/4_0/toolkit/docs/NVIDIA_CUDA_ProgrammingGuide_4.0.pdf
Flynn, M. J., 1972, “Some Computer Organizations and Their Effectiveness,” IEEE Trans. Comput., 100, pp. 948–960. [CrossRef]
Tasora, A., and Anitescu, M., 2011, “A Matrix-Free Cone Complementarity Approach for Solving Large-Scale, Nonsmooth, Rigid Body Dynamics,” Comput. Methods Appl. Mech. Eng., 200, pp. 439–453. [CrossRef]
Taylor, R. L., 1999, FEAP, a Finite Element Analysis Program: Version 7.1 j User Manual, Department of Civil and Environmental Engineering, University of California, Berkeley.
Hoberock, J., and Bell, N., 2010, “Thrust: A Parallel Template Library,” Available at http://thrust.googlecode.com
Gerstmayr, J., and Irschik, H., 2008, “On the Correct Representation of Bending and Axial Deformation in the Absolute Nodal Coordinate Formulation With an Elastic Line Approach,” J. Sound Vib., 318, pp. 461–487. [CrossRef]


Grahic Jump Location
Fig. 1

Upper part: beam with three elements characterized by the generalized coordinates e1 through e4. In the approach adopted herein, the elements in a beam are separated; each element has a set of 12 generalized coordinates associated with it (lower part), which leads to the presence of some simple kinematic constraints (Eq. (11)). This decouples the problem and increases opportunities for parallelism at the price of a larger albeit sparser problem.

Grahic Jump Location
Fig. 2

Advancing the numerical solution in time draws on three loops: the outer one marches forward in time with step size h; at each time step tn+1, the second loop solves the algebraic discretization problem in Eq. (18) using a stopping test controlled by the parameter εNS. Each iteration in this second loop launches a third loop whose role is that of producing the correction δn+1(k) in the value of the acceleration and Lagrange multipliers. The corrections are computed at various degrees of accuracy by controlling the value of the residual εNS(k) in the stopping test of the BICGSTAB iterative solver [17] applied for the linear system in Eq. (22).

Grahic Jump Location
Fig. 3

Displacement in the x-direction of the pendulum tip (ANCF, abaqus, and FEAP comparison)

Grahic Jump Location
Fig. 4

Displacement in the y-direction of the pendulum tip (ANCF, abaqus, and FEAP comparison)

Grahic Jump Location
Fig. 5

Displacement in the z-direction of the pendulum tip (ANCF, abaqus, and FEAP comparison)

Grahic Jump Location
Fig. 6

Conservation of total energy for the time evolution of a flexible pendulum using the Newmark integration algorithm

Grahic Jump Location
Fig. 7

The y-position of the tip of the cantilever beam subject to a concentrated force of 60 N down the negative y-direction plotted against the number of elements used to discretize the beam. Equilibrium is reached by the numerical solution (continuous line) at −0.1016 m and by the analytical result (broken line) at −0.1019 m. Closer agreement can be obtained by redefining the curvature in Eq. (8) according to the work performed in Ref. [26].

Grahic Jump Location
Fig. 8

The number of iterations in the nonlinear Newton solver is shown as a function of time. The solver typically converges within 2–4 iterations.

Grahic Jump Location
Fig. 9

A system containing 101,025 beam elements constrained together with 640,146 constraints to form a net configuration. The net was constrained along one of the edges and subjected to a gravitational acceleration of 9.81 m/s2 in the negative y direction. An animation of the simulation can be found at http://sbel.wisc.edu/Animations/#47

Grahic Jump Location
Fig. 10

The computational time for net models of several sizes were recorded using several different GPUs. The time it took to simulate ten time steps of the simulation was plotted as a function of the number of beams required to create the net. The scaling analysis went up to 0.5 × 106 ANCF elements.




Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In