Modeling and Simulation of Contacting Flexible Bodies in Multibody Systems

Full text

(1)Linköping Studies in Science and Technology Thesis No. 989. Modeling and Simulation of Contacting Flexible Bodies in Multibody Systems by. Iakov Nakhimovski. Submitted to the School of Engineering at Linköping University in partial fulfilment of the requirements for the degree of Licentiate of Engineering Department of Computer and Information Science Linköpings universitet SE-581 83 Linköping, Sweden Linköping 2002.

(2)

(3) Modeling and Simulation of Contacting Flexible Bodies in Multibody Systems by Iakov Nakhimovski Dec 2002 ISBN 91-7373-559-0 Linköpings Studies in Science and Technology Thesis No. 989 ISSN 0280-7971 LiU-Tek-Lic-2002:62 ABSTRACT This thesis summarizes the equations, algorithms and design decisions necessary for dynamic simulation of flexible bodies with moving contacts. The assumed general shape function approach is also presented. The approach is expected to be computationally less expensive than FEM approaches and easier to use than other reduction techniques. Additionally, the described technique enables studies of the residual stress release during grinding of flexible bodies. The overall software system design for a flexible multi-body simulation system BEAST is presented and the specifics of the flexible modeling is specially addressed. An industrial application example is also described in the thesis. The application presents some results from a case where the developed system is used for simulation of flexible ring grinding with material removal. This work has been supported by SKF Sweden AB, KK-stiftelsens företagsforskarskola i Linköping and the ECSEL (Excellence Centre for Computer Science and Systems Engineering in Linköping) Graduate School.. Department of Computer and Information Science Linköpings universitet SE-581 83 Linköping, Sweden.

(4)

(5) Acknowledgments Many people contributed to the work presented in this thesis and essentially made the thesis possible at all. First, I would like to thank my supervisor Dag Fritzson at SKF Sweden AB who has suggested many of the ideas realized in this work, encouraged me to tackle the problems during the whole period of studies, and gave many important comments on the text of the thesis. Fredrik Berntsson from MAI, Linköping University deserves a special thanks for helping me in understanding of some of the mathematical problems discussed in the thesis. I really appreciate the reviews of the thesis done by Peter Fritzson, the head of PELAB (Programming Environments Laboratory) at Linköping University, and by Claus Führer, associate professor in Numerical Analysis at Lund University. Theirs comments helped to significantly improve the text. Many thanks go to my fellow industrial PhD students at SKF Alexander Siemers and Aleksandar Bogdanovski for sharing the exciting experience of studying and working on BEAST at the same time. I also want to thank all the other members of the BEAST team at SKF: Mikael Holgersson, H˚akan B˚astedt, Lars-Erik Stacke, Leonid Gershuni, and Alexei Jolkin as well as Daniel Lindgren from Xdin for many fruitful discussions and lots of useful feedback. Special thanks to Lars-Erik Stacke for proof reading the thesis. Many thanks also to all my colleagues at PELAB who have created a nice working atmosphere at the University. A special thanks goes to Bodil Mattsson-Kihlström for handling all the administrative work caused by my frequent travels to Gothenburg. Furthermore, I am grateful to all the administrative staff at IDA, in particular Lillemor Wallgren, Britt-Inger Karlsson, and Bodil Carlsson. The work was financially supported by the SKF Sweden AB, ECSEL (Excellence Center for Computer Science and Systems Engineering in Linköping) Graduate School and KK-stiftelsen företagsforskarskola in Linköping. Finally, I would like to thank my parents in Russia, and my relatives and friends in different countries around the world for constant moral support and belief in my ability to do the work and write this thesis.. Iakov Nakhimovski Linköping, November 2002.

(6)

(7) Contents Notation 1. 2. 3. 4. Introduction 1.1 Simulation of Mechanical Systems 1.2 The BEAST Toolbox . . . . . . . 1.3 Model Requirements . . . . . . . 1.4 Related Work . . . . . . . . . . . 1.5 Thesis Overview . . . . . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. 8 8 9 12 14 15. Key Equations 2.1 Flexible Body Motion and Shape Functions . . . . . . . . . . . . . 2.2 Generalized Newton-Euler Equation . . . . . . . . . . . . . . . . . 2.3 Mass Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Quadratic Velocity Vector . . . . . . . . . . . . . . . . . . . . . . . 2.5 Generalized Forces . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.1 Generalized Elastic Forces . . . . . . . . . . . . . . . . . . 2.5.2 Generalized Viscosity Forces . . . . . . . . . . . . . . . . . 2.5.3 Generalized External Forces . . . . . . . . . . . . . . . . . 2.5.4 External Point Force . . . . . . . . . . . . . . . . . . . . . 2.5.5 External Point Moment . . . . . . . . . . . . . . . . . . . . 2.5.6 External Volume Load . . . . . . . . . . . . . . . . . . . . 2.5.7 Generalized Forces from the Internal Stress Release during Grinding . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5.8 An Interpolation Method for Forces with Discontinuities . . 2.6 Rotation of a Material Particle due to Deformation . . . . . . . . . 2.7 Using Intermediate Coordinate System . . . . . . . . . . . . . . . . 2.8 Calculating Jacobian . . . . . . . . . . . . . . . . . . . . . . . . .. 17 17 19 20 22 23 23 24 27 27 28 29. General Shape Functions 3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Choice of Shape Functions . . . . . . . . . . . . . . . . . . . . . . 3.3 Separation of Elastic and Rigid Body Motion Modes . . . . . . . . 3.4 Special Shape Functions for Solid Bodies in Cylindrical Coordinates. 40 40 41 42 44. 1. 30 31 34 34 36.

(8) 3.5 3.6 3.7. Imposing Constraints for Rigid Body Motion . . . . . . . . . . . . Reducing the Number of Flexible States . . . . . . . . . . . . . . . Volume Integration . . . . . . . . . . . . . . . . . . . . . . . . . .. 45 46 48. 4 Modeling Ideal Connections with Control Points 4.1 Control Points and Flexible Bodies . . . . . . . . . . . . . . . . . . 4.2 Boundary and Loading Conditions . . . . . . . . . . . . . . . . . . 4.3 An Example of Control Points Usage . . . . . . . . . . . . . . . . .. 51 51 53 54. 5 Static, Eigenmode and Quasi-static Single Body Analysis 5.1 Static Loading Cases . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Eigenmode Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Quasi-static Analysis . . . . . . . . . . . . . . . . . . . . . . . . .. 56 57 57 58. 6 Verification 6.1 Long Beams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 A Thin Ring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 60 60 61. 7 Simulation System Design 7.1 System Overview . . . . . . . . . . . . . . . . . . . . . 7.2 Rigid Body Model Classes . . . . . . . . . . . . . . . . 7.3 Flexible System Design . . . . . . . . . . . . . . . . . . 7.3.1 Limitations of the Rigid Body Model Design . . 7.3.2 Implementation of the Control Point Architecture 7.3.3 Object-Oriented View of Shape Functions . . . . 7.3.4 Flexible Body Class Design . . . . . . . . . . . 7.4 A Simple Model Example . . . . . . . . . . . . . . . . 7.5 Solution Procedure . . . . . . . . . . . . . . . . . . . . 7.5.1 RHS Evaluation . . . . . . . . . . . . . . . . . 7.5.2 Jacobian Evaluation . . . . . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. . . . . . . . . . . .. 63 63 64 66 66 66 68 68 69 70 70 72. 8 Dynamic Simulation of Grinding 8.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . 8.2 Need for Grinding Simulation . . . . . . . . . . . . . 8.3 Grinding Machine Model . . . . . . . . . . . . . . . . 8.4 Modeling Tribological Contacts with Material Removal 8.5 Grinding simulation results . . . . . . . . . . . . . . . 8.6 Summary . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. . . . . . .. 73 73 73 74 75 76 81. . . . . . .. 9 Conclusions. 82. 10 Future Work. 84. A Shapes Generated by General Functions. 85. 2.

(9) B Example Acceleration Calculations for a Flexible Ring Bibliography. 88 103. 3.

(10) Notation The notation used throughout the thesis is described below. The section can be used as a reference when reading the equations presented in the thesis. The notation used here is both the direct matrix notation of vectors and tensors, and the vector-dyadic notation. Therefore, the basic quantities and some formulas are given in both notations. The orthonormal right-handed rectangular Cartesian coordinate system is used. For further reading we refer to [9, 7, 11, 8] or other well-known standard publications.. Variables & basic definitions. . . transformation matrix from coordinate system 1 to system 2.. . a scalar.. . rules can be applied..

(11) a vector or first-order tensor on which vector transformation. a unit vector.. . . a skew symmetric matrix constructed from the components of 3-element vector . That is a matrix:. . The important properties of a skew symmetric matrix are:. " !. . T. # $ ( T. normally a second-order tensor or matrix, in some cases an array.. &%'

(12) %. T. a second-order tensor or matrix.. transpose of. . . and . 4.

(13) . (c). ) * (,+ +.-0//. a 3-element vector expressed in component form in coordinate system ’c’.. /. c. time derivative of a vector with respect to coordinate system c. c. If no coordinate system is specified then the global inertial coordinate system is assumed.. ! 1 1 ! 1 !87. :

(14) ; . 2 %435%. (or. '6 ! 1 ) dot product of two vectors. cross product of two vectors. 2 %43

(15) %' 2 %43

(16) %59, 9. dot product of two matrices. a constant that gives the relative stiffness between object 1 and object 2.. 7 83 <'=>3 -? A ( A ( AB ( ADC. 3 -? . damping matrix.. ff. a control point. defined in the coordinate system. unit base vector in a rectangular Cartesian axis system Euler parameters. E. E , GE F H

(17) ;. elastic (Young’s) modulus.. J J , K M9 L (OJ N 9PL Q. an identity matrix. I. matrices of elastic and viscous material constants. a external force refering to object 1. ext. R. moment of inertia for object 2 relative to coordinate system 1 some of the inertia shape integrals defined by the equations 2-9, 2-10 and 2-13 respectively. stiffness matrix.. ff. S. 3@< .. length or width. S

(18) ; S. a moment relative to the origin of system c refering to object 1. 1/c. an external force couple, i.e. moment, refering to object 1. ext. a force couple, i.e. moment, refering to object 1 and caused by interaction with object 2. 1,2. S SUTVT ( X S WYW ( S. complete symmetric mass matrix of a body as used in NewtonEuler equations. ff. S. sub-blocks of complete mass matrix . The subscipts specify the degrees of freedom (DOF) corresponding to the particular sub-block. ’ ’ specify the three translational DOFs, ’ ’ the three rotational DOFs and ’ ’ elastic DOFs. The nondiagonal. Z. [. 5. \^]. A.

(19) S_T`W. blocks of the matrix use mixed subscripts, e.g., stands for the sub-block defining inertia coupling between the translational and rotational DOFs.. ab (Mc e. a/b. e fe. (b) f (b) f. d. a “point”, i.e. the vector between a and b. If b is a coordinate system, then it refers to the origin of system b. a ”point” vector expressed in the coordinate system b.. g9. fl T ( fl ( fl g ng m g l W lp o l rq s . i. wx9 . j. k. total generalized force tensor as appears on the right hand side of the Newton-Euler equation defined in 2.2. f. parts of a generalized force vector corresponding to translational, rotational and elastic DOFs respectively. external generalized force tensor. quadratic velocity tensor. generalized force due to internal residual stress release. the location of the origin of coordinate system 2 relative to the origin of coordinate system 1. u\^] t 9. 0t v tN t N 9PL. . h. an element of the point vector associated with the axis . Axes are named ’ ’, ’ ’, ’ ’ or ’1’, ’2’, ’3’ in different equations as appropriate in the particular context.. l. t. d. mass and mass density of object . The subscript is omitted if the discussion clearly indicates a single object.. a shape matrix, that is a 3 matrix representing the deformation field of a flexible body. is used to denote the -th row of the matrix, that is the row associated with the -th Cartesian coordinate axis.. h. h. a dynamic inertia shape vector defined by the Equation 2-15. one of the inertia shape integrals defined by Equation 2-11. nine inertia shape integrals defined by Equation 2-12. The subscripts and specify corresponding coordinate axes.. h. ?. e. f < (zy{( f i ( j ( k g f}| (z~( g k g. components of the deflection vector f(b) in the direction of coordinate axis . Where is one of or for a Cartesian coordinate system and one of for a cylindrical coordinate system.. h. f. h. the tensor of the generalized elastic coordinates of a flexible body. volume. velocity. 6.

(20) ( ( W * a. angular acceleration of ’a’ relative system 1 differenciated in system 2. If no system 2 is specified then system 1 is assumed.. . (M a. . potential energy, strain energy and work of external forces.. . rotation angle. a first order tensor containing three angles representing the rotation of coordinate system 2 relative to system 1. This tensor cannot be expressed in a coordinate system. the angular speed of coordinate system 2 relative to coordinate system 1 . the angular speed of “a” relative to coordinate system 1.. d -th vibration mode frequency [radians/second].. elastic strain tensor. viscose stress tensor. elastic stress tensor. Lame’s constants. Poisson’s ratio. infinitesimal rotation angles.. 7.

(21) Chapter 1 Introduction 1.1 Simulation of Mechanical Systems Simulation technology becomes more and more important to industry. Simulation systems can be seen as virtual test rigs that allow for replacement of long and costly experiments on a real test rig with faster, cheaper, and more detailed simulations on a computer. Mechanical systems are obvious candidates for simulation for many reasons, some of them are listed here:. . Building a prototype mechanical system just to test a new design is very expensive and time consuming. Sometimes manufacturing of a small slightly modified part at a factory or workshop requires several weeks.. . Many machines may be dangerous to use without simulation in advance. Simulation nowadays is a must, e.g. an airplane or car design.. . Some measurements in a real dynamic system are very hard or even impossible to perform. However in a simulation all the variables are accessible.. . Tuning the system parameters is much easier for a simulated system than in the real machine. Some particular effects that in real systems are obscured by some other phenomena can be isolated and carefully studied in a simulation. On the other hand, in other situations such small effects can be completely neglected to allow more careful study of the main relations.. Historically, two main types of simulation of mechanical systems have been dominating:. . Multi-body dynamics simulations typically deals with systems of interconnected rigid bodies. In many cases the connections are limited to idealized 8.

(22) . joints (e.g., revolute, prismatic, etc) and traditional force elements (e.g., spring, damper, external load). Multiple dynamic contacts are often important parts of the model and a relatively simple contact model is most often employed for performance reasons. Multi-body systems are mostly used for time domain dynamic simulations. Finite element tools are mostly used for static and quasi-static (with constant angular velocities) analysis of systems and components as well as for modal analysis of structures. They are also used for detailed simulations of very small and specific parts of a system, e.g., a single contact. Finite element analysis generally requires more computational effort than multi-body analysis.. The general development of simulation software goes in the direction of simulation of complete systems with more and more detailed and accurate analysis. This means that multi-body systems now include deformable bodies and more detailed analysis in the models. At the same time finite element tools are beginning to be used for the dynamic simulations. This thesis contributes to this process on the multi-body systems side. It discusses the floating-frame of reference formulation which is gradually becoming the standard way to simulate flexible bodies in the context of a multi-body simulation. To further narrow the field of the research we limit our development to multi-body systems specializing in detailed contact analysis. An example of such a system is the BEAST (BEAring Simulation Tool) described in the next section. The toolbox is used industrially at SKF - the world’s leading rolling bearing manufacturer. An industrial simulation system is always the result of interdisciplinary research where contributions from different fields are integrated into a single generally useful system. In the case of flexible multi-body simulation system contributions from mechanics, numerical analysis, and computer science are necessary. Here, mechanics provide the mathematical models of physical phenomena in the form of equations; numerical analysis comes up with the methods and algorithms for the solution of the equations, and computer science gives the design guidelines for the efficient implementation of the model on a computer. The aim of this thesis is to describe an approach for building a flexible multibody simulation system as an extension to the existing rigid body one (BEAST). Therefore all aspects of the system development have to be considered. The thesis discusses the necessary mathematical model for deformable bodies simulation, provides some important numerical algorithms, and discusses the object-oriented design of the computer implementation.. 1.2 The BEAST Toolbox BEAST (BEAring Simulation Tool) is a fully three-dimensional simulation program developed at SKF and originally used to perform simulations of bearing dynamics 9.

(23) on any bearing type. This has enabled studies of internal motions and forces in a bearing under any given loading condition. The bearing could be loaded in any way the user required. BEAST has been seen as a virtual test rig where the user has full insight into the dynamic behavior of the bearing components. The toolbox development in terms of design generalization has widened the application area of the program. Later versions of BEAST include models of grinding machines, experimental engines, transmissions, etc. BEAST can now be seen as a general multi-body simulation tool specialized in detailed contact analysis.. Beauty Input Output. Input Animations Output. Input. Input. ViewBeast. BEAST. RunBeast Output. 2D plots Output. Output. Input. Out2In. Figure 1-1: Programs in the BEAST toolbox. The BEAST toolbox includes a set of tools that are interacting with each other by means of different kinds of files. Figure 1-1 shows the programs in the toolbox and their interactions with the file storage. The tools are: Beauty is an advanced 3D graphical tool designed for setting up the model and simulation parameters in the input files and visualizing the output files. Some of the features of the tool include animating the simulated sequence from different viewpoints, visualization of force vectors and surface associated data (such as pressure distribution). The visual representation of the simulation results in Beauty in an easily understandable way contributes to a quicker interpretation of the result and popularity of the complete toolbox among the users. Figure 1-2 shows a snapshot of the Beauty working on an input file for a grinding machine model. ViewBeast is specially designed for 2D plot presentation and analysis. In addition 10.

(24) Figure 1-2: A snapshot of the Beauty displaying a grinding machine model. to the basic functionality, i.e., curve plotting, different operations on the simulation results are possible. For example, the user can apply Fourier transforms to a variable or specify his/her own function on several variables. RunBeast is a remote simulation interface system. Its architecture is described in [22]. The tool provides a user-friendly interface for submitting a simulation on a remote computation server which is often a parallel computer. BEAST is the main simulation program. It reads in a model specification from input files, performs the simulation and generates a set of output files containing the results. Most of the developments of this thesis were realized within the framework of this tool. Out2In is a small utility that allows generation of new input files from the simulation output files. In this way the user is given a chance to interrupt the simulation, modify some parameters, and continue the simulation from the moment where the previous run was terminated. This means that is the results of one 11.

(25) simulation can be used as the initial conditions for another. From the mathematical point of view the BEAST system solves the simultaneous system of Newton-Euler equations of motion for every body in the mechanical system. The Newton-Euler equations are formulated as second order ordinary differential equations (ODEs) on explicit form. The second order differential equations system is rewritten as a first order system. Typical characteristics of such ODEs are: mathematical stiffness, very high numerical precision of the solution required by the application, and computationally expensive evaluation of the derivatives. See Section 2.8 and [13] for more details. The tool uses a modified version of the CVODE [4] differential equation solver for the numerical solution of the resulting system of equations. This solver is one of the most well known free implementations of the implicit backward differentiation formulas (BDF) integration scheme. The most important features of this solver are:. . A variable order multi-step method that requires high continuity of the derivatives.. . Adaptive time step changes based on a local error estimate. Efficient use of the solver requires a fast procedure for the system Jacobian calculation (i.e., the matrix of partial derivatives of the state derivatives with respect to the state variables).. The original solver has been extended to handle special kinds of Jacobians (block diagonal with borders and sparse) and to perform a simultaneous RHS and Jacobian calculation [13]. For the work described in this thesis, the BEAST system has been the assumed implementation platform. All the ideas for the modeling, the algorithms, and the design decisions were tried for applicability and usability in this system. Hence we can say that the system has become a test implementation of the presented approaches. Its successful industrial usage validates the correctness of the decisions taken.. 1.3 Model Requirements The main aim of this thesis is to provide the elastic body model suitable for dynamic simulations involving multiple moving contacts. The main issue is an efficient model of the overall structural deformation. Local deformation in a contact is covered by the contact model discussed in other reports on the BEAST (BEAring Simulation Tool) and GRIT(GRInding simulation Tool) systems [6, 23, 12]. The requirements on the model can be subdivided into three equally important groups. They are: the computational performance or efficiency, the accuracy, the potential to incorporate the analysis needed for grinding. 12.

(26) . . The computational complexity of a flexible system model depends mostly on the number of state variables used for the flexible body. Hence, an efficient model should use as few states to describe flexibility as possible. The system of ODE (ordinary differential equations) that is generated in the dynamic contact analysis problems is in most cases mathematically stiff. The commonly used type of numerical integrator for such systems is multi-step variable time step solver. Important consideration for such solvers is the length of the time step in the integration. Longer time steps means less computations and shorter simulation times. That is why an efficient flexible body model for the case of moving contact should not force the numerical solver to take significantly smaller time steps compared to the same model with rigid components. That is, the higher frequencies in the system are expected to be associated with the contact forces calculations and not with the flexible body eigenfrequencies. Let us provide a very simple example to illustrate the need for accuracy. In Figure 1-3 a rotating elastic ring pressed between two plates is shown. From the physical point of view the radial stiffness of the ring is independent of its orientation. Hence assuming the constant force acting on the ring the elastic deflection, measured as the distance between the plates, should become constant after some transient process. When simulating this in the time domain, it is important that we do not introduce vibrations due to modeling errors or numerical errors. The situation can get worse for a more complex model where small, induced vibrations can occasionally resonate with some resonance frequency of the system leading to completely unrealistic results.. . Figure 1-3: A rotating flexible ring pressed between two plates. In the third group of requirements one can mention distortion, due to the initial residual stresses in the surface layer combined with uneven material removal. 13.

(27) 1.4 Related Work A 3-D FEM model can be used to model flexible bodies. The required procedures are described in a wide variety of articles and books. References [1, 5] provide extensive discussions of the methods. However, finite element methods for transient analysis often suffer from long (days, weeks) computation times. Also, the radial stiffness computed at the element nodes will slightly differ from the stiffness between the nodes. This is due to the fact that the deflection between finite element nodes is normally computed using a low order interpolation scheme. For our simple example described above, variation in radial stiffness leads to induced non-physical vibration with the frequency dependent on the number of elements on the circumference of the ring and rotational speed. Without doubt the effect can be minimized and made negligible by using a large number of finite elements but that would lead to an even longer computational time. The effects of the discretisation errors in dynamic simulations are discussed in [16]. A commonly used approach to efficiently simulate structural elasticity is to employ some kind of reduction technique [5, 29, 30] from a detailed finite element model. Reduction schemes often dramatically decrease the number of state variables used in dynamic simulation. The reduction is done prior to the simulation using a FEM-tool. The result of a reduction procedure is a set of assumed displacement shape functions that is used in the dynamic simulation. The main challenge of a reduction with modal synthesis is to select a small set of superimposed mode shapes that give acceptable calculation precision for the simulation. It is one of the most complicated and discussed problems in flexible multi-body simulation (see [20]). The mode shapes are calculated using data from static mode shapes and/or frequency response modes, and eigenmode shapes. These resulting mode shapes can be calculated only for some specific boundary conditions, i.e., given attachment or interface nodes. However, in the case of a moving contact the boundary conditions are changing dynamically. Interaction between any two points of the contacting surfaces is possible. Hence many nodes on the surfaces have to be marked as interface nodes. Large number of interface nodes leads to the large number of state variables and makes the dynamic simulation slower. The reduced model still has the same vibration problem as the FEM model, due to the discrete interface nodes on the contact surfaces. The solution to this is to employ continuous shape functions for the interface motion, i.e., control the surface interface nodes by the shape functions. This can be done directly in the reduction scheme [30], or may be done as a separate transformation step after a standard reduction with all interface nodes. These reduced models cannot handle the requirements from the third group mentioned above in Section 1.3, i.e., non-linear effects due to initial stresses and material removal. 14.

(28) 1.5 Thesis Overview The dynamic equations of motion for the rigid bodies in a multi-body system can be defined in terms of the mass of the body, the inertia tensor, and the generalized forces acting on the body. On the other hand, the dynamic formulation of the system of equations of motion for linear structural systems requires the definition of the system mass and stiffness matrices as well as the vector of generalized forces. This thesis summarizes the equation of motion for deformable bodies that undergo large translational and rotational displacements using the floating frame of reference formulation. The set of the inertia shape integrals required for the equations is defined. These inertia shape integrals, that depend on the assumed displacement field, appear in the non-linear terms that represent the inertia coupling between the reference motion and the elastic deformation of the body. Some numerical procedures required to perform a dynamic simulation involving flexible bodies including numerical volume integration and Jacobian matrix evaluation are also described. There are several kinds of analysis that are possible only for flexible body models. This thesis describes three general types of such analysis: free-body eigenmodes analysis, body deformation under static load and under quasi-static conditions. Results of these kinds of analysis can be used to quickly assert the basic properties of the flexible body model before running computationally heavy dynamic simulation. There is also a discussion on a more specialized kind of analysis: ring distortion due to the residual stress release during grinding. The thesis also provides some examples that were used to verify the correctness of the mathematical model presented. The main intention of this thesis is to provide a complete and possibly compact description of the equation and procedures necessary to perform a dynamic simulation of a flexible multi-body system and some simpler kinds of single body analysis. The thesis is therefore system implementation oriented and does not generally cover the mathematical derivation of the presented equations. The complete derivation of these equation can be found in most books on flexible dynamics (see [21]). The thesis is structured into seven chapters and two appendixes.. . . . The first chapter provides an introduction and an overview of the thesis. The second chapter introduces the mathematical model of an elastic body using the floating frame of reference formulation and assumed shape function. All the important equations that define the model are listed and explained in this chapter. Appendix B shows how the equations presented in this chapter can be used for modeling of a flexible ring. The third chapter presents the set of general shape functions that are proposed to be used as assumed mode shapes for elastic bodies of used in simulations with detailed contact analysis. The necessary requirements and boundary conditions as well as volume integration procedure applicable to the presented 15.

(29) . shape functions are discussed in this chapter. Appendix A provides a complimentary example for this chapter.. . Chapter four covers some design decisions that were made to provide an efficient user interface to for specifying boundary and external loading conditions in simulations involving flexible bodies.. . The fifth chapter describes some extra types of analysis that are possible for elastic body models.. . Chapter six gives some simple example test cases that were used to verify the accuracy of the flexible body modeling using floating frame of reference formulation and general shape functions.. . . Chapter seven focuses on the system architecture and design for a multi-body simulation system that includes flexible components. We specially address the problem of transition from a multi-body system that allows only rigid bodies to the system with both kind of components. Chapter eight describes an application example for the simulation system discussed in the thesis. The importance of the structural flexibility and its influence on the simulation results for the particular application are discussed. The last two chapters summarize the results of the performed work and discuss the possible directions for the continuation of the research and development in this area.. 16.

(30) Chapter 2 Key Equations This chapter introduces the mathematical model of an elastic body using the floating frame of reference formulation and assumed shape function. Appendix B shows how the equations presented here can be used for modeling of a flexible ring. The chapter is strongly based on the results of [21] and therefore no derivation for the most of the presented equations is given. The relations that are developed in this thesis and cannot be found in the reference above are:. . Representation in the different coordinate systems.. . Generalized viscosity forces and damping modeling.. . Generalized external moment.. . Generalized external body load.. . Generalized force from the residual stress release. An interpolation method for forces with discontinuities.. . Jacobian calculations. Explicit step-by-step calculation procedure as presented in Appendix B.. In the following sections the subscript that normally indicates the body number will be omitted with the understanding the all vectors and matrices are associated with some particular single body. Since the equations involve only two coordinate systems - the global and body ones - one corresponding letter (’g’ or ’b’) will be used to specify the coordinate system where the vector is expressed or differentiated.. 2.1 Flexible Body Motion and Shape Functions A rigid body in space has six degrees of freedom that describe the location and orientation of the body with respect to the fixed frame of reference. On the other hand, 17.

(31) structural components such as beams, plates, and shells have an infinite number of degrees of freedom that describe the displacement of each point on the component. Since it is computationally impossible to deal with infinite spaces, classical approximation methods are employed wherein the displacement of each point due to elastic deformation is expressed in terms of a finite number of coordinates:. L9 fe g fe 9 g O fe g b¡9 (b) f (b) f (b) f. e. 9 [ 9 , where [ 9 [ 9 f e g £ ¤¤¥ 2 9D'9 , where 9 9 f e g ¤ 3z9B¢.9 , where ¢.9 ¢.9 f e ¦ ¤ g (b) 0. (b) 0. (2-1). (b) 0. e e. where f(b) gives the displacement of an arbitrary point that has coordinates 0(b) in the undeformed state. The vector of displacement (or deformation vector) f(b) is space- and time- dependent. The coefficients , , and are assumed to depend only on time. The above equations can be written in the following matrix form:. 9 2 9. t. [ 9 ( '9. e. ¢¨9. (b) f. \§] `9 2 9. t ! . 3z9. (2-2). f. where is the three-rows -columns shape matrix whose elements are the functions , and ; and f is the vector of elastic coordinates that contains the time dependent coefficients , , and . The total number of elastic coordinates (which is of course equal to the number of columns of the shape matrix ) should be determined heuristically depending on the simulation accuracy required and the importance of the elasticity effects for the complete model. By using the outlined approximations, the global position of an arbitrary point of the body can be defined as. 3z9. \§]. ©. e. (g) P/g. s. (g) b/g. p«P¬ ! e . ª. (b) P/b. (2-3). where. s ®«M¬ e. (g) b/g. (b) P/b. define the origin of the body reference; orthogonal rotation matrix;. e ª e. e ª t ! . (b) (b) is equal to 0(b) f 0 in the body coordinate system. ©. f. and gives the displacement of the point. Note that it is generally possible to specify the shape matrix in different coordinate systems and transform it using rotation matrix and normal rules for coordinate transformation. However, use of the body reference frame enables the most natural and efficient choice of the shape functions. For this reason the shape matrix is always assumed to be calculated in the body coordinate system. Hence the coordinate system needs to be specified for all the tensors in the following sections. 18.

(32) 2.2 Generalized Newton-Euler Equation The generalized Newton-Euler equation for the unconstrained motion of the deformable body that undergoes large reference displacement is given by:. S_TVT p«P¬ t 6v ®«M¬ t N s ± JON W ] ¯ aa°A J0| N WYW3 S ± - d j «M³ f l W ² «M³ l´ ² ¬µ³ f l W g ² ¬µ³ l´ ² f l W Q ! g m 7 ! ) l m g. (g). b/g (b). ff. . Where. . . R. f. l" ² «M³ µ¬ ³ l´l m ² . f. ff. f. ff. (2-4). f l o ² «M³ f l o g ¬µ³ ª f l o g m² g . ff. R. (2-5). f. Alternatively using body coordinate system for all vectors:. S_TVT t 6v t N s ± l" ² ¬µ³ ² ¬µ³ O J } W W J O W " l N N ¯ aa°A | 3 S ] ± l m - d j. . (b). ff. . . f. R. f. R. . b/g (b). R. f. f. . (2-6). where. s ± * b/g. . is the acceleration of the offset of the body reference coordinates relative to. g. the global coordinate system. (b). ± «M¬. f. is the angular acceleration vector of the reference of the deformable body defined in the body coordinate system.. . second derivative of the vector of time-dependent elastic generalized coordinates of the body f . - orthogonal rotation matrix, which for the case of Euler parameters is defined as:. < y f A y f AB y f A A ABADC y f A AB A ADC «M¬ y f A A g ABzADC g < y f A y f ABg y f A AB ª A ADC gg ( g g g ª y f A AB A ADC g y f A AB ª A ADC g < y f A g y f A g (2-7)

(33) ¶&¶ the four Euler rotation parameters, defined for the rotation about where A C are # vector ¸ · ¸ · ( ¸ · ( ¸ · $ on angle ~ as A ¸ · ¹Mº¼»G½ ( A ¸ · ^¹Mº¾»¿½ ( (2-8) AB ¸ · §¹5º¼»G½ ( ADC ÁÀDÂ ¹Ã½ Ä 19.

(34) The other parameters of the equation are the mass matrix, the generalized forces and quadratic velocity vector. The components of the mass matrix are defined in the next section, generalized forces are discussed in Section 2.5, and the quadratic velocity vector that includes the effects of Coriolis and centrifugal forces is defined in Section 2.4. The Newton-Euler equation given in this section can be further simplified by using the shape functions that satisfy mean-axis condition. See Section 3.3 for further development and an efficient solution approach.. 2.3 Mass Matrix Even though the mass matrix is highly nonlinear, the deformable body inertia can be defined in terms of a set of constant inertia integrals that depend on the assumed displacement field.. J ÆÅÈÇ c e ¬µ³ ¬µ³ 9K PL ÅÈÇÊcOË Ì² ; 9 Ë Ì² ; L É t N ÁÅ ÇÊc t t N 9PL Å@Ç c t 69 t L É µ ¬ ³ J N 9PL ÆÅ Ç cÍË Ì² ; 9 t L É. e. (b) 0. É. . ( h ( ? < (y{( É ( h ( ? < z( y{( ( h ( ? < (zy{( . (2-9) (2-10) (2-11) (2-12) (2-13). where (b) 0. . c. is the undeformed position of an arbitrary point on the deformable body; is the mass density of the body; is the volume of the body;. t 9. t. h. is the -th row in the body shape matrix , defined by Equation 2-2.. Note that in the special case of rigid body analysis the shape integrals are given by Equations 2-9 and 2-10 only. The mass matrix components now can be defined: -. S_T{T ÏÎ Ç c I a I É. I. ¿ . , where - is the identity matrix for the spatial case, m - mass of the body. For the case of conservation of mass this matrix is the same for both cases of rigid and deformable bodies. 20.

(35) -. t v. Ð v ; Ð v ; ÐÐ vv ;; Ð v ; Ð v ; t0 v # Ð v ; ( Ð v ; ( Ð v ; M$ J t N ! ] ª. is given by. where -. JON W}W. ¿ symmetric matrix with¬µ³the elements ¬µ³ ¬Y; given by¬µ³ ¬Y; ¬µ³ ¬Y; f Ë Ñ² µ¬ ³ ¬Y; f Ë Ñ² ¬µ³ ¬Y; Y ¬ ; Ë Ñ² Ë Ñ² Ë Ñ² Ë Ñ² g ª g µ ¬ ³ µ ¬ ³ f Ë Ñ² ¬Y; f Ë Ñ² ¬Y; Ë Ñ² ¬µ³ ¬Y; Ë Ñ² ¬µ³ ¬Y; ÅÇ c g ª g ¯ aa°A | 3 f Ë Ñ² ¬µ³ ¬Y; f Ë Ñ² ¬µ³ ¬Y; - d j g ª g is a. (2-14). (2-15). . É. . (2-16). Let us now show how this matrix can be calculated efficiently by using precalculated inertia shape integrals. First we will analyze the components of the diagonal elements:. ¬µ³ Î Ç c r# f Ë Ñ² ¬Y; 9 g $ É Î Ç ÎÇ ÎÇ. ¬µ³ c r# f Ë Ì² ; 9 ª t 9 ! ] g $ É ¬µ³ ¬µ³ c #rf Ë Ì² ; 9 g ª 6] ! t 69 ! t 9 ! ] ª y ! Ë Ì² ; 9 ! t 9 ! ] $ É ¬µ³ c f Ë Ì² ; 9 g É ª 6] !@Î Ç c # t 69 ! t 9z$ É ! ] Ç c # Ë Ì² ¬µ; 9³ ! t 94$ ! y Î ] ª É K 9P9 ª yGJ N 9P9 ! ] ª 6] ! t N 9P9 ! ]. (2-17) Using the procedure similar to the one shown above we can come up with the following compact expressions for the diagonal and non-diagonal members of :. JON W}W. Ò. K Ò W}W ; 9P9 L K LÓL ª y f L J N LÓL g ! ] ª K W}W ; 9PL f K LÓ9 ª f J N LÓ9 ª J N 9PL g ! ] where h ( ? < (zy{( VÕ ?× Ö h. 6] ! f L t N LÔL g ! ] t ª 6] ! N LÓ9 ! ] g. (2-18). JN WYW. Note that for the case of rigid body analysis the elastic coordinates are zero and so the matrix is constant. -. JON W ]. has three rows defined as follows. 6] ! f t N t N f J N J N JON W ] 6] ! f t N t N g f J N J N g 6] ! f t N Y t N

(36) gg ª f J N Y J N

(37) gg 21. (2-19).

(38) -. S ]z]. is independent on the generalized coordinates of the body and, therefore, constant: (2-20). S ]z] t N n t N n t N n ª ª. 2.4 Quadratic Velocity Vector The quadratic velocity vector can be defined as. l o #rf l o 6T f l o 6 f l o 6] $ 6 g gm g. In the three-dimensional analysis the components of the vector. pl o. (2-21) are defined as. £ ¤¤¥ f lp o T² «M³ b«M¬ ! #Øf t0N v y t N ) $ g g ª f lp o f JN WYW ! J)N WYW ! f JN W ] ! ) ¤ (2-22) g g ¤ gm f l o ] Î Ç c.Ù t 6 #Øf ! e ! e ) $}Ú ¦ y ! Á g g ª É the angular velocity vector defined in the body coordinate system. and where is a skew issymmetric matrix given by (2-23) f. (b). (b). (b). (b). f. (b) f. (b). (b). (b). (b) 3. (b) 3 (b) 2. (b) 1. (b) 2 (b) 1. The quadratic velocity vector includes the effect of Coriolis and centrifugal forces as nonlinear functions of the system generalized coordinates and velocities. in terms of the shape integrals we need to evaluate the In order to express symmetric matrix:. f . f lp o ] g. . g fMf . f g g fMf f f g ª g g g fMf f g ª g ª g g (2-24) o f l The g ] part in the explicit form is given below: Ò f l o ] Û Û f % f J N %6 t N &% ! y ! Û Û % Ð % ) (2-25) g g ª r %5 g r %M Ò r %5 % Ð &% ) The term can be further simplified to f t N Y t N

(39) g ! ) ª f t N t N g ! ) ª f t N t N g ! ) (2-26) (b) 3. (b) 1 (b) 1. (b) 2 (b) 3. (b) 2. (b) 1. (b) 1. (b) 2. (b) 3. (b) 2. (b) 3. (b) 1. f. (b) 1 (b) 3. (b) 3 (b) 2. f. f. (b) 3. f. (b) 1. f. 22. (b) 2. f. (b) 2.

(40) 2.5 Generalized Forces 2.5.1 Generalized Elastic Forces The virtual work due to elastic forces can be written as. . . Ü Å Ç × 6 Ü É. (2-27). Ü Ý. is the virtual work of the where and are the stress and strain tensors, and elastic forces. Since the rigid body motion corresponds to the case of constant strains and since we defined the deformation with respect to the body references, there is no loss of generality in writing the strain displacement relations in the following form: where Þ. ÏÞ e . (b) f. (2-28). is a differential operator defined according to:. ß % y< f ¨w ; % ª `w % ; ª Û xw 9 ; ¼wx9 ; % g ( w¨ ; % + w¨% ( d (à < (y{( 9 +i f}| or, in case of cylindrical coordinate system with coordinates (z~( k g ß n + w¨| á + ß n |< f + w ~ w¨á + ½ ª g w ß n + +k ß Y |< f + w¨~ á w g ª + w | ½ + + ß + w¨á + w | ½ +k ª + ß |< + w ~ + w + ª + k½. (2-29) (see [24]):. (2-30). For a linear homogeneous isotropic material, the constitutive equations relating the stress and strains can be written as. E # ß n ( ß n ( ß n ( ß Y ( ß ( ß P$ 6 , &# â n ( â n ( â n ( â Y ( â ( â P$ 6 , E where . (2-31) is the. symmetric matrix of elastic coefficients, introducing the material properties. Substituting Equation 2-31 into Equation 2-27 we can get:. Ü Ý ! Q T f. where. Q. ff. ff. ! Ü . f. ÏÅ Ç f Þ t g 6 E Þ t É 23. (2-32) (2-33).

(41) Q. ff is the symmetric positive semidefinite stiffness matrix associated with the elastic states of the body. Boundary conditions discussed in Section 3.3 ensure that the stiffness matrix becomes a positive definite one. Since the virtual work of the elastic forces does not depend on the rigid body states the corresponding components of the generalized force vector are zero. The resulting vector is therefore given by:. l´ ² «M³ ¬µ³ l´² l m Q ! R. ff. f. f. To complete the description we show how the matrix of Lame’s constants:. E. (2-34). can be expressed in terms. . ª y y ª E ª y 8y 8y y . . The Lame’s constants are related to the engineering constants ulus and - Poisson’s ration with the following equations:. ä¤ f y å¤ ã ª ª g æ y f ª g. . . ã. ä¤ å¤ y f ã ª < g æ f < y ã f < g ª g. (2-35). - Young’s mod-. (2-36). It is important to mention that the stresses and strains we were describing above are the second Piola-Kirchhoff stress and Green-Lagrange strain tensors. The most important observation about these two measures is that the components of the tensors do not change under rigid body translation and rotation and so they provide a natural material description for the case of large rigid body motion but small elastic deflection and strain analysis. The observation is of practical importance since Hook’s law is applicable only to small strains and because in the particular problems we are interested in the condition of large rigid body motion, accompanied by small strains holds.. 2.5.2 Generalized Viscosity Forces We start the description by mentioning a big similarities between the elastic and viscosity forces and the related laws. For instance, for a linear homogeneous isotropic material, the constitutive equations relating the rate of strain change and viscosity stress can be written as (2-37). EGF ) 24.

(42) ç ”Mass proportional” damping. ç ¡. ”Stiffness proportional” damping. ¡. Figure 2-1: Damping coefficient as a function of frequency.. Q. 7. And so following the procedure outlined in the previous subsection we can compute the damping matrix ff similarly to the stiffness matrix ff . However, the material is often not so well known for damping and so we will describe constants’ matrix one of the approximation techniques - Rayleigh damping. This damping model is discussed in more details for the case of finite-element analysis in [1]. The damping matrix in this case is calculated as:. EF. 7 Áèé! S. è. ê. ff. ff. Q ªëê !. (2-38). ff. where (1/sec) and (sec) are constants to be determined from two given damping ratios that correspond to two unequal frequencies of vibration. Assuming that for two and we know the damping ratios and respectively vibration frequencies (in fractions of critical damping), the Rayleigh coefficients can be found from the equations:. . . ç. ì è è ªª . ! ê y ! ç !ê y !ç. ç. (2-39). Similar to the previous section the resulting generalized force vector can be formulated as:. lé ² «M³ µ¬ ³ l"² l m R. Note that given the è. and. ê. . . . . . ) 7 ! ff. f. (2-40). f. coefficients one can establish the damping ratio that 25.

(43) is specified at a value of. í :. ! . ç è î ª y ê í. (2-41). The general view of this curve is shown in Figure 2-1. Simple analysis of this curve tells us that damping coefficient has its minimum value at the frequency and that the coefficient grows linearly for higher and hyperbolically for lower frequencies. In fact, one of the features of Rayleigh damping is that the higher modes are considerably more damped than the lower modes, for which Rayleigh constants have been selected. The solution of the linear system Equation 2-39 gives the following relations for the and coefficients:. ç ðï é è !ê ¡. è. ê. òñ ôè ó ê ¡. ä¤ è y ! f ç ! ç ! g å¤ æ y fç ! ç ! g ê . è. (2-42). When choosing the damping ratios and frequencies it is important to make sure that the resulting and coefficients are positive. The negative value of one of the coefficients leads to the negative damping ratio for a range of flexible body frequencies. Thus the system will become non-dissipative and generate vibrations during the simulation. Since all the parameters used ( , , and ) are positive the following conditions follow from Equation 2-42 if we assume :. ê. ç ç. O õ å äæ ç ! ç ! 0õ ç ! ç ! 0õ ç ®ö ç ÷ö . (2-43). or, alternatively:. (2-44). In practice, the values of the coefficients are often selected based on the experience with simulation of similar structures when the results of a simulation were fitted to the same experimental results. That is the same and are used in the analysis of similar structures. The magnitude of the Rayleigh coefficients is to large extent determined by the energy dissipation characteristics of the construction, including the material. Another use of Rayleigh damping might be to improve the numerical performance of the simulation by damping high frequency vibrations outside the frequency band that the user is interested in. For example, suppose that only one damping ratio for the first body eigen-frequency is known. Then, it might be reasonable to select some high frequency that limits the frequency band and set the corresponding damping ration to one thus damping out all the high frequencies. Having the two damping ratios and for two distinct frequencies the Equation 2-42 can be used to obtain the and coefficients.. è. çW. è. çø çW. ê. ùø çø. W. 26. ê.

(44) 2.5.3 Generalized External Forces The virtual work of all external forces acting on a body can be written in compact form as (2-45). Ü éW l 6W Ü . l W. where is the vector of generalized external forces associated with the body generalized coordinates. In a partitioned form, the virtual work can be written as. Ü s Ü éW #rf l 6T f l W 6W f l W 6] $ Ü ú (2-46) g g g Ü f l T and f l W W are the generalized forces associated, respectively, with the where g g f l W ] is translational and rotational coordinates of the selected body reference, and g (g). (g) e. f. (g) e. the vector of generalized forces associated with the elastic generalized coordinates of the body. Generalized forces may depend on the system’s generalized coordinates , velocities, and on time. Later in this section in order to provide simpler and more compact equations the small angles and actual moments vector are used instead of generalized rotational coordinates (Euler parameters) and associated part of the generalized force . That is we will use:. . f l W 6W g. . fl W ú gm. Ü s Üé W #rf l 6T f l W 6 f l W 6] $ Ü g g m g Ü (g). (g) e. (2-47). f. 2.5.4 External Point Force. H. (g). H f ( - g (g). Let us assume that a force acts at point body. The virtual work of this force is defined as. e. (g). of the deformable. Ü éW f H 6 Ü e (2-48) g where e is the global position of point © of the deformable body. Let the position . of the point be defined by the vector e e s ª b«M¬ ! e (2-49) is the local position of point © with respect to the body coordinate system. where e Ü is then defined as The virtual change e # «P¬ ! e $ ! Ü «M¬ ! t ! Ü Ü e Ü s ª + + ª (2-50) Ü s «M¬ ! e ! Ü «P¬ ! t ! Ü ª (g). (g) P/g. (g) P/g. (g). (g) b/g. (g) P/g. (b) P/b. (b) P/b. (g) P/g. (g) P/g. (g) b/g. (b) P/b. (g). f. (b) P/b. 27. f.

(45) The vector. Ü e. (g) P/g. Ü e. can be written in a partitioned form as. Ü s # I ( P« ¬ ! e ( «M¬ ! t $ Ü Ü (g). (b) P/b. (g) P/g. where the shape matrix as. t. (2-51). f. defined at point. ©. . Thus the virtual work. Ü éW. is defined. Ü s Ü éW f H 6 # I ( M« ¬ ! e ( f M« ¬ ! t $ Ü ú g Ü g (g). (g). (b) P/b. (2-52). f. and so the generalized forces in Equation 2-47 can be recognized as. T H ( f l W e ! µ¬ }« ! H ( f l W ] t 6 ! ¬µ}« ! H (2-53) g ng m g eî e 6 «M¬ µ¬ }« were used to simplify the equations. where the facts that and 6 fl. (g) e. (g). (g). (b) P/b. 2.5.5 External Point Moment. (g). S S. (g). f ( g. e ©. Let us assume that an external moment is applied at point of the deformable body. Let the position of the point be defined by the vector (g) . The virtual work of this moment is defined as. Ü éW f S . Ü . (g). g6. Ü . (g). (2-54). (g). is a virtual change of the the global orientation of the point where vector rules apply to the infinitesimal rotations one can write:. Ü . Ü . (g). Ü ª (g) r. «M¬ ! Ü . (b) f. ©. . Since. (2-55). Ü . (g). (b). where r is the virtual change of the reference frame orientation and f is the virtual change of the material particle orientation inside the flexible body as defined in Section 2.6. In terms of the generalized coordinates the virtual orientation change is then defined as:. Ü . (g). Ü s # ( b«M¬ ( p«M¬ + $ Ü + Ü (b) f. (g). f. f. 28. (2-56).

(46) û û. (b). pÝ\§]. where f is a matrix calculated by differentiation of Equation 2-68 with f respect to the elastic states:. + t + k + < + t y +i + +t +j (b) f f. t #t 6 ( t 6 ( t 6 $6 . + t + t j + +k +t +i. (2-57). ©. is the shape matrix defined at point . Finally the resulting generalized force can be recognized as:. fl. T fl g ( W ng m «M¬ where the fact that 6 (g) e. ¬µ}« ! S ¬µ}«. (g). ( f l W g ] f + g 6 ! ¬µ}« ! S + (b) f. (g). (2-58). f. was used to simplify the equations.. 2.5.6 External Volume Load In this subsection we will analyze the external loads acting on all of the particles of the body. The most common cases of such loads are gravity, magnetic and electrostatic forces. For generality we will assume that some load vector with the . dimension N/m (Newton per cubic meter) is acting on each volume particle For the case of gravity this load is given as , where is the gravity acceleration constant vector and is the mass density of the body. Using the equations developed in Section 2.5.4 and integrating over the volume we can get the components of the generalized body load :. ü. . c ! ý. c. ý. É. . lÿ þ f l° þ T ÈÅ Ç ü ² «M³ ü ² «P³ ÅÈÇ ü ² «M³ ! (2-59) g É É ² µ¬ ³ ! b¬µ}« ! ü ² «P³ Å Ç e ² ¬µ³ ! ¬µ}« ! ü ² «M³ < t v ! ¬µ}« ! ü ² «M³ f l° þ Å Ç e c gnm É É (2-60) f l° þ ] ÅÈÇ t 6 ! ¬µ}« ! ü ² «M³ Å@Ç t 6 ! ¬µ}« ! ü ² «M³ < t N ! ¬µ}« ! ü ² «M³ c g É É (2-61) t0 v and t N were defined in Section 2.3. For the case of the gravity force where the equations can be further simplified considering that a c ! : f l T a ! ýô ² «P³ f l g t0 v ! ¬µ}« ! ý «P² «M³ ³ (2-62) f l gm] t N ! ¬µ}« ! ýô ² g 29.

(47) Note that the derivation as presented in this section can easily be generalized for other kinds of loads that do not depend on the elastic states. That is as long as the load (force) vector does not depend on the point position within the body one can sum (or integrate) over the points once at the start up and use the precalculated shape integrals during the simulation. One useful application of this approach is area load, where a constant load acts over a surface area. From the user interface design point of view, the body and area loads rise an issue of sub-volume and sub-area specification in the input. That is a simulation tool needs to provide some means to define the geometry to be integrated. Both for the area and volume specification two approaches are possible. Closed boundary A sub-volume can be defined by a set of surfaces and a sub-area - by a set of curves on the surfaces. In both cases the boundary must be closed to make integration possible. Parametric specification The sub-volume and sub-area can be specified by coordinate or parametric intervals. This approach is essentially a special case of the closed boundary specification where the boundaries have simple configuration. For some cases, e.g., when the specified sub-volume is completely inside the body, the parametric specification simplifies the integration procedures since the integration ranges become constant.. 2.5.7 Generalized Forces from the Internal Stress Release during Grinding The structural elasticity model presented in this thesis was incorporated into a grinding simulation tool [12]. One of the special effects that is of interest for this simulation is the internal stress release during grinding of rings. The internal residual stresses in the material can be caused by material phase changes in the hardening process. When stressed material is removed the released stresses result in elastic distortion of the body. These distortions are especially significant for large rings where they often cause manufacturing problems and therefore should be incorporated into the flexible model [10]. The internal stresses do not affect the elastic behavior of the body as long as the material is intact. When the material is removed they result in additional generalized force that affects the generalized elastic coordinates but does not influence the rigid body part of the system. Following the procedure described in Section 2.5.1 we can get: (2-63). l q ÁÅ Ç × 6 Þ t « ( É. . «. where the volume integration is done over the material volume that is ground away during the simulation and stands for the residual stress tensor. For simulation purposes it is often possible to assume that all the surface particles are subjected to the same treatment and only a thin layer of material is removed 30.

(48) . during grinding. That is why it is realistic to model the internal stress tensor as a constant tensor in the surface layer when the local coordinate system (s) normal to the surface is used to calculate the strains and stresses. Assuming that a rotation matrix transforms the body coordinates into the local surface coordinates one can perform integration of the shape matrix independently:. ¬. ³ l ² q

(49) f ² ³ 6 Å Ç Þ ² ³ f ¬Mt « ´ g gÉ ¬. (2-64). f | (~( k. g. Note that for the case of ring grinding where a cylindrical coordinate system is used as the body coordinate system the matrix does not depend on the angle . The problem with the evaluation of the generalized force comes from the fact that the volume representing the ground material changes with time. In a simulation the ground surface is described by a mesh. Every node of this mesh ground at the point on the surface. The keeps the value of the height integration in Equation 2-64 should obviously be performed over this surface mesh. During a grinding simulation the values are updated only once for every ODE solver time step and are kept constant when the solver does trial time steps during internal solver iterations. It means that the generalized force can also be evaluated only once per time step and can not be accurately evaluated between the are accumulated during the simulation run and can not be steps. The values directly computed from the state variables. Therefore the extrapolation procedures inside the solver can not predict the change of the resulting force. Hence, from the ODE solver point of view the force changes discontinuously. For a high-order solvers with adaptive time step size like CVODE force discontinuity results in convergence problems and shorter time steps, leading to much longer simulation times. In order to avoid this discontinuity problem an interpolation procedure was developed. The procedure is actually general and can be used for smoothing of any kind of force input that causes numerical problems. The procedure is discussed in Section 2.5.8. ~. l q. «. ¢ « fw (¸ g. ¢ « fw (¸ g. ¢ « fw (¸ g. fw (¸ g. l q. 2.5.8 An Interpolation Method for Forces with Discontinuities Some non-linear effects modeled in dynamic simulations make it impossible, or very hard, to evaluate the force function at every time instant. In such cases some extrapolation or interpolation procedure must be used to provide a continuous force function in time domain. High order ODE solvers used for dynamic simulations assume high continuity of all the functions involved in the problem. Some experiments done in the BEAST tool show that at least second time derivative’s continuity is needed in order to get reasonably large time steps in the solution. The interpolation procedure described here can be used for the forces which are changing with lower frequencies than some other fast processes in the simulated system. The main idea of this procedure is to delay the real force value for a small time 31.

(50)

(51) f- g -. -. Interpolation between and. f-g. - ª Üv. f- , ) f- ,± f- g g g f Figure 2-2: Force interpolation as af filter block. A continuous force function - g generated for the discreet input

(52) - g .. . is. interval and to use cubic interpolation to get a smooth change of the force function between the current and the delayed values. Before using the procedure the interpolation time delay parameter must be selected. The only requirement on this parameter is to be larger than the maximum possible time step of the ODE solver. Initial value of the force variable to be interpolated must be given. The initial first and second time derivatives of the variable are normally assumed to be zero if no other information is available. The interpolation procedure can be seen as a filter block shown in Figure 2-2. The only input for the block is the calculated force . This value is delayed and which is seen as the expected force value for time . assigned to the variable So, at every time , such that the interpolation block is using the following four variables:. Üv. -. ]. - - ö f - ª Ü v g.

(53) f- g. - ª Üv. f- g; First and second time derivatives ) ) f - g and ± ± f - g ; Delayed value of force ] . . Force variable current approximating value. f-g. Having the four parameters above (two values and two derivatives) it is possible to calculate the approximate value by using cubic interpolation:. f - g f < Ü < v f - )f f ª - - g <. - g g ª Ü < v f -. ] Ü < v f - - g - ± y <Ü v f - - f - - Ü v g g g g 32. (2-65).

(54) . ] . . . ¬ . ] -. -. -. - - ¬ - ª Ü v - ¬ . - ª Ü v. . Figure 2-3: One dimensional force interpolation. The values of the force variable for the time instances and are calculated using third order interpolation between and for and between and for .. . ]. -. . ] . . f-g.

(55) f- g. Note that is a continuous function which is the essential output of the interpolation block and should be used instead of the discreet values . When the numeric ODE solver has done some iteration and decided the step the force variable value can be size and state variables values for the time calculated. In order to move to the next time step and to interpolate the values for times between and we need to evaluate and the derivatives and according to the equations:. - . ) . - f - g ± f f ) - Ü v - g ª ] - Ü v - g ) f < f - Ü - g ± f f < f - Ü v - g -g v y ª ª f f ± - Ü v - g ª ] - Ü v - g ) f - Ü v - g ± f < f - Ü v - g g ª - ] - ) . ] . (2-66). Figure 2-3 shows an interpolation scenario where the ODE solver evaluates one extra time instance between two time step and . The corresponding value of the force variable is interpolated from , , and using the procedure above.. -. . 33. ± .

(56) The interpolation procedure described in this section provides continuity of the second time derivative of the force vector which is sufficient for most cases. It is certainly possible to achieve higher order continuity by using higher order interpolation methods. The interpolation procedure can be applied element-wise to vector variables and hence can be used for generalized forces’ interpolation.. 2.6 Rotation of a Material Particle due to Deformation In general during the deformation of a body, any element is changed in shape, translated, and rotated [24]. In order to analyze systems with rotational springs and dampers it is necessary to be able to calculate the orientation of a material point in the deformed body. In order to calculate the orientation of a material particle for the case of large deformations one has to use the properties of the deformation gradient [1]. However for the case of small deformations where only linear terms are taken into account the rotation can be analyzed as infinitesimal. The rotation matrix that performs the rotation of a particle from the undeformed state to the infinitesimally deformed orientation can be calculated as (see [21]):. ] Ì I (2-67) ª f I where is an identity matrix and ( ( g are three small angles corresponding to the rotations around the three respective coordinate axes. The angles can be calculated as: + e + e +j +k < + e + e (2-68) y +i + k + e + e +j +i (b) f,y. (b) f,z. (b) f,z. (b) f,x. (b) f,x. (b) f,y. 2.7 Using Intermediate Coordinate System It is sometimes convenient, and/or more efficient, for the numerical solver to observe the motion of a body in a coordinate system other than the inertial one. For example if the simulated body is a part of a complex system and it is desirable to analyze the body movements from a certain known position in the system. From the mathematical point of view, this known position is a coordinate system with known motion. In such a case a special care must be taken when working with time derivatives of 34.

(57) position and orientation of a body. Below, the transformation rules for taking time derivatives with respect to different coordinate systems, are given. To be more specific let us assume that two coordinate systems are involved in the analysis.. . . System is the inertial coordinate system where the Newton-Euler equations are valid.. s! ". #s ) " " *. $s ± " "*. System is a coordinate system with know motion in . That is, the position vector , linear speed , and linear acceleration as well. % ) " " *. ". as rotation matrix acceleration. Then for a physical vector. ". , relative angular velocity. and relative angular. are known.. s. s ) " * s ) * " s ª. the following transformation is valid: (2-69). Applying Equation 2-69 on the relative angular velocity vector shows that the time derivative is invariant with respect to these two systems:. ) " " * ) " * " " ) " * (2-70) ª Therefore the coordinate system where differentiation ) " is done can be omitted from . the specification the angular acceleration vector Using the same Equation 2-69 twice we arrive to the relation for the second time derivative:. s ± " * s ± * % ) " s y % " s ) * % " f % " s (2-71) g ª ª ª Now we will consider a body coordinate system 3 defined relative to the inter . That is the position of the body is defined with the mediate coordinate s and itssystem vector orientation with a rotation matrix . Hence the position of.

(58). & is defined as s

(59) " s$ " s

(60) . the body in the coordinate system.

(61). ª. (2-72). Then Equation 2-69 leads to the following expression for the linear velocity of the body: (2-73). s )

(62) " " * ! s ) " " * s )

No results found