Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

(1)

UPTEC F 19032

Examensarbete 30 hp Juni 2019

Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

Christoffer Fjellstedt

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

Christoffer Fjellstedt

Two methods using lumped element (lumped parameter) methods to model flexible bodies have been presented. The methods are based on the

concept of using a Guyan reduced stiffness matrix to describe the elasticity of a body. The component to be modeled has been divided into two parts using FE software and the mass and inertia tensor for the respective part of the component have been retrieved. The first method has been based on including the elements from the stiffness matrix in compliant constraints. The compliant constraints have been derived and a prototype has been implemented in MATLAB. It has been shown that using compliant constraints and stiffness parameters from a Guyan reduced stiffness matrix it is possible, with highly accurate results, to describe the deformation of a flexible body in multibody simulations.

The second method is based on springs and dampers and has been implemented in the simulation environment Dymola. The springs and dampers have been constructed to include coupling elements from a Guyan reduced stiffness matrix. It has been shown that using the proposed method it is possible, with highly accurate results, to describe the static deformation of a flexible body. Further, using dynamic simulations of a full robot manipulator model, it has been shown that it is possible to use the spring-damper model to capture the deformation of the links of a manipulator in dynamic simulations with large translations and rotations.

Examinator: Tomas Nyberg Ämnesgranskare: Robin Augustine Handledare: Hans A. Andersson

(3)

1 Popul¨ arvetenskaplig sammanfattning

För i princip alla maskiner, anläggningar och liknande är det nödvändigt att utveckla korrekta modeller som beskriver systemen. Dessa modeller kan användas för att styra, kontrollera och simulera ett system eller en maskin. Ett omr˚ade där kontroll och simuleringar är viktigt är inom robotik, där man vill kunna kontrollera robotens rörelser men ocks˚a genom simuleringar kunna analysera en robots beteende med hög noggrannhet. Utvecklingen inom industrirobotar har g˚att mot mer lätta strukturer med fördelar s˚asom lägre produktions- och energikostnad. För lättare strukturer blir emellertid betydelsen av elasticiteten av de olika komponenterna av roboten mer betydande, vilket ställer krav p˚a att elasticiteten inkluderas i de modeller som beskriver roboten.

Det finns flera möjliga sätt att modellera elasticiteten p˚a en komponent. Ett sätt är att använda s˚a kallade finita element (FE) modeller. Denna typ av modell beskriver en komponent genom att dela upp strukturen i ett ändligt antal sm˚a strukturer. Idén är att det är enklare att beräkna dessa sm˚a komponenter var för sig och sedan sätta samman detta till ett fullständigt system än att försöka beräkna den fullständiga komponenten direkt. FE modeller kan ofta beskriva en komponent eller ett system med hög noggrannhet men är vanligen beräkningsmässigt för krävande, dvs. det tar l˚ang tid att lösa systemen.

En annan möjlig väg för att inkludera elasticitet i modeller är att använda en s˚a kallad lumpade element (lumpade parameter ) metod. Detta sätt att modellera elastiska komponenter g˚ar ut p˚a att man placerar fjädrar i axlarna mellan stela kroppar. Axlarna kan vara faktiska komponenter p˚a en robot eller fiktiva, dvs. introducerade endast för att möjliggöra simulering av elastiska komponenter. En elastisk kropp, t.ex. en balk, kan delas upp i flera stela kroppar och elasticiteten introduceras i systemet via fjädrar i de fiktiva axlarna. Lumpade element metoder resulterar ofta i snabbare beräkningstider än fullständiga FE modeller.

I denna uppsats har tv˚a metoder för att modellera elastiska komponenter med hjälp av lumpade element presenterats. Den första metoden bygger p˚a att begränsa rörligheten mellan tv˚a komponenter och att i denna begränsning införa en elasticitet (compliant constraints). Den andra metoden bygger p˚a att direkt använda matematiska uttryck för fjädrar och dämpare för att introducera elasticitet i systemet.

För b˚ada metoder har samma generella modelleringskoncepts använts. Den elastiska komponent man har avsett att modellera har delats upp i tv˚a delar med hjälp av s˚a kallad FE mjukvara. Fr˚an samma mjukvara har även styvhetsparametrar för komponenterna exporterats, dvs. värdet p˚a fjädrarna. Dessa parametrar har för en komponent skapats med en modellreduceringsprocess som kallas Guyan reducering.

Denna metod bygger p˚a att man reducerar ner en fullständig FE modell till en betydligt enklare modell som beskrivs av färre frihetsgrader. Guyan reducering är en statisk reduceringsmetod, vilket betyder att den kommer att beskriva den statiska deformationen p˚a en komponent korrekt men inte den dynamiska.

Det har visats i uppsatsen att b˚ada metoder med hög noggrannhet kan beskriva den statiska deformationen av en enkel balkgeometri, en mer komplex geometri samt ett mer komplext system av tre sammansatta balkar. Metoden baserad p˚a fjädrar och dämpare har ocks˚a testats som delar av en fullständig robotmodell i dynamiska simuleringar, där det har visats att modellen har god potential att användas i simuleringar av fullständiga robotmodeller.

(4)

2 Introduction

In order to predict and control the behavior of a machine, industrial robot, vehicle and any other device it is necessary to develop correct and accurate models of the systems. The trend in robotics has been towards more light weight structures, where the aim is to increase the performance of the robot while decreasing production costs [23, 21]. A consequence of the move towards more light weight structures is that the systems become more sensitive to the elasticity of the links, which have effects on both the static deformation and the vibration modes of the robot. To enable correct motion control of a robot it can therefore be necessary to include the elasticity of the links in the models. However, if the model is going to be used in control or simulations it is necessary that the complexity of the model is kept to a minimum. Complex models are often capable of producing very accurate descriptions of the robot systems, but this comes at the cost of significant computation times. The aim is therefore both fast and accurate simulations of complex elastic behavior.

AGX Dynamics is a physics engine developed and maintained by the Ume˚a based company Algo- ryx Simulation [3]. The physics engine is based around rigid body systems with joints and non-linear dynamics (collisions, contacts with friction and so on). In this modeling approach simple and complex systems, e.g. robots can be simulated in real time with high accuracy. Elasticity can be introduced in the systems by using the lumped element method, in which the elasticity of a body is represented by splitting a body into multiple rigid bodies and inserting joints with 6 degrees of freedom between the bodies. The elasticity of the joints can for a simple geometry, e.g. a beam, be calculated from linear elasticity theory. ABB Robotics use AGX Dynamics for internal research and development and in the simulation and programming environment RobotStudio.

The challenge consists in that the actual geometries of robot components can not easily be estimated with simple beam geometries. It is therefore often necessary to use more complex models to describe the different components of a robot. This can be achieved using the finite element method (FEM), which supports dynamic simulations. The resulting system is however too complex, with extensive computation times, to be used directly in dynamic simulations. Therefore, the complex FE models are often reduced using a model order reduction process. This approach is currently used by ABB Robotics in dynamic simulations of flexible link models [38, 37]. The reduced models are however still too complex to enable very fast simulations, or at least simulations near real-time.

The aim of this thesis is therefore to investigate concepts on how the stiffness parameters from reduced FE models can be introduced in lumped element models. This is studied in two parts: firstly, a modeling approach based on compliant constraints is developed, which can potentially be introduced in AGX Dynamics. Secondly, a model based on spring-damper pairs is implemented in the existing simulation environments at ABB Robotics.

2.1 Goals of the thesis

The goal of the thesis is to investigate methods and develop concepts on how to use stiffness parameters from reduced FE models to simulate flexible components in multibody simulations. The goals can be broken down in:

• Study the general concepts of model order reduction of FE models and the possibility to use the reduced models in lumped element methods.

• Develop a general model for simulations of flexible structures using a lumped element approach.

• Develop a prototype of a model based on compliant constraints in MATLAB, which potentially can be implemented in AGX Dynamics. Evaluate the performance of the implemented prototype model on a simple beam geometry and a more complex geometry

• Develop a prototype of the modeling concept in the simulation environment Dymola. Evaluate the performance of the implemented prototype model on a simple beam geometry, a more complex geometry and in dynamic simulations with a full robot manipulator model.

The thesis is limited to mainly study quasi-static deformations. This means that the models will be developed with the aim to describe the deformation of structures for primarily static situations with the possibility for accurate results for low frequency situations.

(7)

2.2 About the project

This thesis has been part of a cooperation project between ABB Robotics in V¨aster˚as, Sweden, and Algoryx Simulation in Ume˚a, Sweden. The parts of the thesis concerning the general modeling concepts and the solution based on compliant constraints (section 4.4) has been the main results of the collabora- tion. The parts of the thesis concerning the spring-damper model implemented in Dymola (section 4.6) has solely been performed at ABB Robotics.

ABB Robotics is a world leading manufacturer of industrial robots and has operations in 53 coun- tries. Over 300 000 ABB robots are installed over the world. ABB Robotics supply complete robot products, with systems and services. Key markets include automotive, metal fabrication, foundry, solar and consumer electronics. [1]

Algoryx Simulation is a leading provider of software and services for physics based interactive simulations. Algoryx was started in 2007 as a spin-off from Ume˚a University in Sweden. Algoryx provides products such as AGX Dynamics and Algoryx Momentum. [3]

(8)

3 Theory

In this section the theory behind fundamental concepts used in the project will be presented.

3.1 Modelling of robot manipulators

In this section an overview of basic concepts and the state of the art regarding modeling of robot manipulators will be given. In section 3.1.1 basic concepts will be presented. In the following sections an overview of the state of the art of dynamic models of robot manipulators will be given.

3.1.1 Basic concepts

The definition of what exactly constitutes an industrial robot is somewhat debated. Often what is meant with the term industrial robot is a device that can perform a wide variety of tasks. A distinction is in that case made to fixed automation devices which mostly are limited to performing a single task. [6] There exist several types of structures that are commonly called industrial robots. The two main typologies are serial chain and parallel mechanisms. The serial chain system is a linkage of bodies, where all bodies, except for the first and last member, is connected to two other bodies. A parallel structure consists of two or more parallel chains of bodies, where all chains connect to the same platform or end effector. [30]

Figure 1 shows the robot manipulator that will be considered in the scope of this thesis. The manipulator consists of seven links connected by six revolute joints that allow relative rotation between neighboring links. The number of degrees of freedom (DoF) of a manipulator are the number of indepen- dent position variables that need to be specified in order to locate all parts of the manipulator. For an industrial manipulator with an open serial structure, each joint is usually defined by one variable, and therefore, each joint equals one DoF. [6] This means that the manipulator in figure 1 is described by 6 DoF. The figure shows the rotational axis of each joint and the commonly used names of the different links. The non-moving base is the first link in the serial chain. After which, the first joint (axis 1) connects the base to the second link, the frame. This continues until the seventh and last link, which is the turning disc.

Figure 1: Industrial robot manipulator. [2]

(9)

At the free end of the serial link, i.e. the turning disc, a tool can be mounted. This part of the manipulator is usually called the end effector. The specific tool depends on the application of the robot and can, for example, be a gripper, wielding torch or electromagnet. [6]

The term tool center point (TCP) is often used to describe the position of a robot. The TCP can be a point on the robot or the end effector. Often the TCP is the focal point of the tool, e.g the tip of a welding torch. If the TCP is described in relation to the tool, it is typically necessary to reprogram the robot if the tool is changed. [12, p. 410]

3.1.2 Dynamic models of robot manipulators

Dynamic models of robot manipulators describe the relationship between the motion of the manipulator and the forces that cause the motion. Dynamic models are used in, for example, simulations or real-time control of robot manipulators. [21]

The most basic model of a manipulator is a collection of rigid bodies, in which all links and joints of the robot are modeled as rigid bodies. In this formulation of a manipulator no information about the elasticity of the components are considered. The rigid bodies are described by their mass, position of center of mass and inertia tensor.

A real flexible object is described by a continuous non-linear system of partial differential equations (PDE) with an infinite number of DoF. The complexity of this kind of description of elastic components usually makes it impracticable in real applications. For including elastic components in multibody systems it is therefore necessary to find a description of an elastic component which is based on a finite number of DoF. To achieve this, three approaches are usually identified [21]:

• Finite element models can be used to acquire a very accurate description of the deflection of a flexible body (see section 3.4). In this formulation a large number of DoF are used to describe a structure. This usually results in very accurate results, but at the cost of more extensive demands on the computational strength of the hardware. Therefore, these models are not widely used in dynamic simulations and control applications.

• Assumed modes models can be used to reduce the partial differential equations describing a structure to a finite number of DoF. Using this formulation the deflection of a structure is described by a finite number of modes. This approach is often used in simulation of flexible structures and is used in dynamic substructuring (see section 3.5.2).

• Lumped parameter models in which the elasticity of a structure is described by a set of discrete and localized springs. Using this approach a structure can be divided into a number of discrete rigid elements connected by non-actuated joints. These joints can be called pseudo-joints since they do not describe a real physical joint, but rather represent a modeling concept. The lumped parameter model is a simple and straightforward approach on how to include elasticity in multibody systems. This approach is however associated with difficulties in determining the values on the spring constants.

The aim of this thesis is to develop on the concept of lumped parameter models, and especially on the incorporation of more complex springs with coupling terms (see section 3.4.1) in the models and on the use of stiffness parameters from reduced finite element models as spring constants (see section 3.5).

For dynamic modeling of flexible robot manipulators, two main approaches can be identified: flexible joint models and flexible link models. [23] These two approaches will be briefly be presented in the following sections.

3.1.3 Flexible joint models

Flexible joint models are lumped parameter models and are usually based on the assumption that the main elasticity of a manipulator occurs in the drive train, i.e. in the gear boxes between motors and links. The rigid links are in this type of model connected by a set of springs and the flexible joints are often modeled with linear torsional springs. It is also common to model the joints with spring-damper pairs. The elasticity in the joints can also be used to capture the elasticity of the links. [38]

In [26] and [22] it was stated that the flexible joint model is not able to accurately describe a modern robot manipulator and a new model called extended flexible joint model was introduced. This model has

(10)

three-dimensional spring-damper pairs in joints one to three and one-dimensional spring-damper pairs in the remaining joints. In [23] it was further shown that the flexible joint model with only torsional elasticity in the gearboxes is not able to accurately describe a modern industrial robot manipulator and that the extended flexible joint model is able to produce significantly better results.

3.1.4 Flexible link models

If the elasticity of the links of the manipulator can not be disregarded or localized to the joints a flexible link model can be used. The elasticity of the links are here described using an assumed modes model. The modeling of flexible links have mainly been considered for especially elastic and light-weight manipulators used in research [23]. In [38] a flexible link model is presented where the links are described by reduced FE models (see section 3.5). The implementation is shown able to describe the behavior of a robot manipulator with high accuracy using only development and datasheet data.

3.2 Multibody systems

In this section an overview is given of important areas in multibody systems. An extensive presentation of the subject can be found in [29].

A multibody system can be understood as a collection of rigid or flexible components, often called bodies. The motion of the bodies in the system can be kinematically constrained with joints, and the bodies are allowed to undergo large translations and rotational displacements.[29] An illustration of the concept is given in figure 2, where two rigid bodies are interconnected with a joint. The joint constrains the motion between the bodies by removing one or more DoF between the bodies. Common types of mechanical joints are revolute, spherical and prismatic joints.

Figure 2: Multibody system with two rigid bodies interconnected with a joint.

3.2.1 Rigid body mechanics

A rigid body is defined as a collection of particles, where the relative distance between any two particles remain absolutely fixed at all times and for all configurations. The motion of a rigid body in space is completely described by six coordinates, where three of these coordinates describe the translation of the body and three of the coordinates define the orientation of the body.

Figure 3 shows a rigid body i in three dimensional space. The global position of any point Pi on the body can be defined as

ri= Ri+ Aiu¯i, (1)

where Ri is the position vector of the body frame of reference expressed in the global (inertial) frame of reference, Ai is the rotation matrix from the frame of reference of the body to the global frame of reference, and ¯ui is the position of the point Pi in terms of the frame of reference of the body. Since the body is rigid the distance between the origin Oi of the frame of reference of the body and the point Pi

on the body is constant. The transformation described by the rotation matrix Aiis not unique. Several different representations of orientations can be found in the literature [29, p. 11-12]. This will be further developed in section 3.2.2.

(11)

Figure 3: Rigid body in three dimensions.

The velocity of point P_i can be found by differentiating equation 1 with respect to time d

dtri= d

dtRi+ d

dt(Aiu¯i). (2)

The product rule applied on the last term gives d

dt(A_iu¯_i) = d

dt(A_i)¯u_i+ A_id

dt(¯u_i) = d

dt(A_i)¯u_i, (3)

where it used that ¯uiis constant since the body is rigid, i.e. _dt^d(¯ui) = 0. Using this and switching to dot notation to indicate differentiation with respect to time, equation 2 can be written as

˙r_i= ˙R_i+ ˙A_iu¯_i. (4)

It can be shown that the time derivative of a rotation matrix is [29, p. 47-48]

A˙i= ω^×_i Ai, (5)

and

A˙_i¯u_i= ω_i× u_i= −u_i× ω_i= −u^×_i ω_i, (6) where ω_iis the angular velocity vector of body i expressed in the global frame of reference and ω^×_i and u^×_i are the skew symmetric matrices of the angular velocity vector and position vector u_i, respectively (the superscript × is used to designate a skew matrix). If the coordinate system is assumed to be expressed in Cartesian coordinates, the skew symmetric matrix for body i is defined as

u^×_i =





0 −u^z_i u^y_i u^z_i 0 −u^x_i

−u^y_i u^x_i 0



, (7)

and similarly for the angular velocity vector

ω^×_i =





0 −ω_i^z ω_i^y ω_i^z 0 −ω^x_i

−ω^y_i ω_i^x 0



. (8)

3.2.2 Representing orientation

There exists several different conventions on how to represent the orientation of a body in 3D space.

These are, among others, rotations matrices, Euler angles, Euler parameters, Rodriquez parameters, direction cosines and Quaternions (see for example [29] or [30]). In the scope of this thesis primarily

(12)

Euler angles and fixed angles will be used to represent orientations [30, p. 13-15]. Therefore an overview of important concepts regarding these conventions will be given below.

If two frames of references with the basis vectors (ê^x_i, ê^y_i, ê^z_i) and (ê^x_j, ê^y_j, ê^z_j), respectively, are considered, the rotation matrix describing the orientation of frame i relative to frame j can in general terms be written as

j iA =





r₁₁ r₁₂ r₁₃ r₂₁ r₂₂ r₂₃ r₃₁ r₃₂ r₃₃



. (9)

Further, a rotation of frame i about the ˆe^z_j axis thorough an angle θ can be written as

A_Z(θ) =





cos θ − sin θ 0 sin θ cos θ 0

0 0 1



. (10)

Similarly, the rotation by an angle θ about the ˆe^y_j axis can be written as

AY(θ) =





cos θ 0 sin θ

0 1 0

− sin θ 0 cos θ



 (11)

and the rotation about the ˆe^x_j axis is

AX(θ) =





1 0 0

0 cos θ − sin θ 0 sin θ cos θ



. (12)

To describe multiple consecutive rotations, the rotation matrices can be multiplied as

k

iA =^k_jA^j_iA, (13)

where ^j_iA represents the rotation of frame i relative frame j and ^k_jA the rotation of frame j relative frame k, which if multiplied in the correct order leads to a description of frame i relative frame k.

Eular angles describe the orientation of a body as a rotation about an axis of a moving frame of reference. Taking frame i as the moving coordinate system and frame j as a fixed frame of reference, one sequence of Euler angles can be described as first a rotation of frame i about ê^z_j by an angle α, then a rotation by an angle β about axis ê^y_i of the rotated frame of reference i, and finally a rotation by an angle γ about the axis ê^x_i of the now twice rotated frame i. This type of representation of rotation is often called Z-Y-X Euler angles. In some literature this representation is also referred to as Tait-Bryan angles. Figure 4 shows the rotation of frame i using Z-Y-X Euler angles. For every rotation a ”prime”

sign is added to the axes of the rotated frame.

Figure 4: Z-Y-X Euler angles.

(13)

If the rotation matrices in equation 10 to 12 are used, it can be shown that the final orientation of frame i relative frame j is described by the following rotation matrix

j

iA_Z⁰_Y⁰_X⁰(α, β, γ) = A_Z(α)A_Y(β)A_X(γ)

=





cαcβ cαsβsγ − sαcγ cαsβcγ + sαsγ sαcβ sαsβsγ + cαcγ sαsβcγ − cαsγ

−sβ cβsγ cβcγ



, (14)

where cα = cos α, sα = sin α, and so on, and the ”primes” on to the subscripts are used to indicate that the rotation is described by Z-Y-X Euler angles.

Using the notation presented in equation 9, the Z-Y-X Euler angles (α, β, γ) can be found from the rotation matrix in equation 14 with the following expressions [6, p. 42-43]

β = atan2(−r31, q

r²₁₁+ r²₂₁), (15)

α = atan2(r₂₁ cβ,r₁₁

cβ), (16)

γ = atan2(r32

cβ,r33

cβ), (17)

where atan2(y, x) is used to compute arctan(y/x). The atan2(y, x) function is used, since it uses the sign of both y and x to identify in which quadrant the resulting angle is situated. This information would be lost if arctan(y/x) was to be used. Therefore, if angles of the full range of 360^◦ are to be calculated it is preferable to use the atan2(y, x) function.

Fixed angles describe the orientation by three rotations about a fixed reference frame. If frame j again is taken as the fixed frame and frame i as the frame to be rotated, a sequence of fixed angle rotations is as follows: frame i is first rotated by an angle γ about the ê^x_j axis of frame j, then about ê^y_j by an angle β and, finally, about ê^z_j by an angle α. This sequence is referred to as X-Y-Z fixed angles and is illustrated in figure 5, where again a ’prime’ sign is added to frame i every time it is rotated. This convention of rotation is also often called roll, pitch and yaw angles, but it should be observed that this name is also often given to other similar types of conventions. The rotation matrix describing X-Y-Z fixed angles is

j

iAXY Z(α, β, γ) = AZ(α)AY(β)AX(γ)

=





cαcβ cαsβsγ − sαcγ cαsβcγ + sαsγ sαcβ sαsβsγ + cαcγ sαsβcγ − cαsγ

−sβ cβsγ cβcγ



, (18)

where the ”primes” on the subscript are dropped to distinguish fixed angles from Euler angles. The rotation matrix describing X-Y-Z fixed angles can be observed to be exactly the same as the one describing Z-Y-X Euler angles. This means that the fixed rotation angles (α, β, γ) can be calculated using equations 15 to 17.

Figure 5: X-Y-Z fixed angles.

(14)

3.3 Constrained multibodies

An effective method of modeling multibody systems is to impose kinematic constraints on the bodies of the system. As presented in section 3.2.1 an unconstrained rigid body is completely defined by 6 DoF. A kinematic constraint imposes some type of limitation on one or more of the DoF of the body, e.g. limiting the body to only move along a straight line or to stay on a surface. Constraint based methods of modeling multibody systems is especially useful when the system to be modeled is governed by complicated and maybe badly understood principles, but the geometrical interpretation of the system is conceptually clear. Constrained multibodies have, among other, been used in simulating beams as lumped elements, wires and granular material [24].

In this section will be given a brief overview of the theory behind constrained multibodies, constraint regularization, a specific time integration scheme called SPOOK and a physics engine based on this scheme, AGX Dynamics.

3.3.1 Multibody dynamics

Let each body in a multibody system of N bodies be described by a state vector (q^T, ˙q^T), where qi = (X^T_i, θ^T_i) is a vector containing the generalized coordinates of each body. The vector Xi is the position vector to the center of mass of body i and θi is the orientation of the body described by Euler angles (see section 3.2.2). The state vector of each body in the system is placed in a partitioned vector q = (q1, q2, ..., qN)^T. In a similar way the velocities of each body are organized as ˙qi= ( ˙X^T_i, ω^T_i), where X˙iis the velocity of X and ωiis the angular velocity vector of body i. The velocities are also partitioned in a vector ˙q = ( ˙q1, ˙q2, ..., ˙qN)^T.

The mass m_i and inertia tensor I_i expressed in the global frame, for respective body, are compiled in a mass matrix

Mi=m_iI_3×3 0

0 I_i

, (19)

where I3×3 is an identity matrix with the dimensions 3 × 3. The system mass matrix for the entire system becomes M = diag(M1, M2, ..., MN).

The Lagrangian of a constrained mechanical system with the notation for the generalized coordinates q and the velocity ˙q as defined above, can be written as [24]

L(q, ˙q, λ, ¯λ) = T (q, ˙q) − U (q) − FR(q, ˙q) + λ^Tg(q, t) + ¯λ^T¯g(q, ˙q, t), (20) where FR(q, ˙q) is the Rayleigh dissipation function, T (q, ˙q) =¹₂q˙^TM ˙q is the kinetic energy, and U (q) is the potential energy. The coordinates of the system q are constrained by holonomic constraints g(q) = 0, with corresponding Jacobian G = ^∂g_∂q and Lagrange multiplier λ. The mathematical interpretation of a position constraint is that it restricts the coordinates of the system to move on a surface, which is a subspace of the original configuration space [28]. The velocities ˙q are constrained by non-holonomic constraints ¯g(q, ˙q, t) with the Lagrange multiplier ¯λ. This brief presentation will be restricted to non- holonomic constraints on the form 0 = ¯g( ˙g, t) ≡ ¯G(q) ˙q. [24]

It is common to use some form constraint regularization scheme to ensure numerical stability [24]. An approach to achieve this is to make the constraints explicitly compliant, which means that the system is allowed to oscillate about the constraint hypersurface g(q) = 0, and then let the system be strongly damped by adding physical damping terms to the Lagrangian formulation. [28] The constraints can be introduced as the limit of infinitely strong potentials and dissipation functions [24]:

U(q) = 1

2g^Tg, (21)

FRγ = 1

2γ¯g^T¯g, (22)

where , γ → 0. It can be shown that by using a Legendre transform, a stiff force can be transformed into a constraint with finite regularization [24]

(15)

U= λ^Tg −

2λ^Tλ, (23)

FRγ(q, ˙q) = ¯λ^T¯g −γ 2

λ¯^T¯λ. (24)

The Euler-Lagrange equation of motion for the corresponding system can be written as a differential- algebraic equation (DAE) [24]

M¨q + ˙M ˙q − G(q)^Tλ − ¯G(q)^T¯λ = f , (25)

λ + g(q) = 0, (26)

γ ¯λ + ¯G(q) ˙q = 0, (27)

where the force vector f = ∇qU (q) − ∇q˙FR(q, ˙q) are the forces on the system and ∇q = _∂q^∂ is the gradient.

It is possible to introduce a physically based dissipation by including a damping parameter τ in the Rayleigh dissipation function FRτ= ₂^τg˙^Tg and using a Legendre transform, which results in [24]˙

F_Rτ = τ ˙λ^Tg −˙ τ 2

λ˙^Tλ.˙ (28)

This results in that equation 26 can be rewritten as [24]

λ + τ ˙λ + g(q) + τ G ˙q = 0. (29)

The compliance and damping factors τ and γ in equation 27 and 29 are expressed as scalars. It should however be observed that these factors are not restricted to being scalars. It is possible to express these factors as matrices. [24] In the scope of this thesis the compliance will be assumed to be in matrix form, resulting in the potential energy

U(q) = 1

2g^T⁻¹g, (30)

where the corresponding constraint force can be written as

f_c= G(q)^Tλ = −G(q)^T⁻¹g(q). (31)

3.3.2 Numerical integration

By using a semi-implicit Euler type discretization and linearising the constraints with a Taylor approx- imation g(q + ∆q) ≈ g(q) + G∆q, it can be shown that the system described by equations 25, 27 and 29 can be formulated as the following discrete stepping scheme [24, 19, 18]:





M −G^T − ¯G^T

G Σ 0

G¯ 0 Σ¯









˙ q_n+1

λ λ¯



=





M ˙q_n+ hf_n

−_h⁴Υg_n+ ΥG ˙q_n 0



, (32)

q_n+1= q_n+ h ˙q_n+1, (33)

where h is the time-integration step size and the regularization and stabilization matrices are Σ = diag 4

h²(1 + 4τ /h)

, (34)

Σ = diag¯ γ h

, (35)

Υ = diag 1 1 + 4τ /h

. (36)

The stepping scheme presented in equation 32 is commonly referred to as SPOOK. The precise regularization and stabilization terms can be derived using a discrete variational principle, which helps to

(16)

construct a time-stepper that preserve the essential physical symmetries. The SPOOK scheme is proved to be linearly stable [18]. The perturbation matrices Σ and ¯Σ are used to handle ill-posed or ill-conditioned problems and the stabilization terms −⁴_hΥg_n+ ΥG ˙q_n are used to oppose constraint violations.[24]

It is also possible to impose restrictions on the constraints and Lagrange multipliers. A non- penetration constraint can be introduced as an unilateral constraint g(q) ≥ 0. The instantaneous (non-smooth) change in velocities at impacts can be model by the Newton impact law as an velocity constraint. Dry friction can be modeled with a Coulomb law condition between the multipliers for normal and tangential contact forces, |λt|≤ µ|λn|. This changes the linear system of equations into a mixed complementarity problem (MCP).

In this work only holonomic constraints will be considered. The stepping scheme in equation 32 and 33 can therefore be rewritten as

M −G^T

G Σ

˙qn+1

λ

=

M ˙qn+ hfn

−_h⁴Υgn+ ΥG ˙qn

, (37)

qn+1= qn+ h ˙qn+1. (38)

It is necessary to choose an appropriate value for the damping parameter τ . If small enough time steps h are used, a possible approach is to identify a physical value for the corresponding system. A safe choice for the value of the damping parameter is however τ ≥ 2h [19]. This choice should result in sufficiently high error decay for the solution and will therefore be used in this thesis. Using τ = 2h and 1/(1 + 4τ /h) ≈ h/4τ it is possible to rewrite equation 34 and 36 as

Σ = 1

2h², (39)

Υ = diag1 8

, (40)

where the compliance now is assumed to be a matrix. It should be emphasized that the strong damping acts in the direction orthogonal to the constraint hypersurface g(q) = 0. All motion in the constraint hypersurface is free and undamped.

3.3.3 Jacobian

Only considering holonomic constraints, the Jacobian G corresponding to a constraint can be found by taking the time derivative of the constraint g(q) and using the chain rule

0 = ˙g(q) = ∂g

∂q dq

dt = G ˙q. (41)

3.3.4 AGX Dynamics

AGX Dynamics is a physics engine built on the concept of constrained multibody simulations. It is developed and maintained by the company Algoryx Simulation AB, located in Ume˚a, Sweden. [3] The physics engine is more precisely oriented towards simulations of rigid multibody systems with contacting and non-smooth dynamics. The numerical time integration in AGX Dynamics is based on the SPOOK stepping scheme. AGX Dynamics is equipped with both direct and iterative MCP solvers. It is also constructed with a hybrid direct-iterative solver utilizing an efficient splitting scheme, which enables simulations with high performance, accuracy and scalability. The physics engine can, among others, be used for simulations of robots, vehicles and granular materials. [27, 5]

3.4 Structural mechanics

The description of the mechanics of structures is an important area. It can, for example, be used to predict the bending behavior of beams or to construct complex models of mechanical components in so called FE software. In this section it will therefore first be given a brief overview of the equation of motion for a linear system with n-DoF. Then the finite element method will be presented, and finally the analytical expressions describing the deformation of a beam will be given.

(17)

3.4.1 Linear n -DoF systems

A spring is an elastic component that relates force to deformation. The idealized spring is assumed to be massless and it can for most structural materials be assumed to follow a linear relationship for small deformations

fs= ke, (42)

where k is the spring constant describing the stiffness of the spring and e = u2− u1, i.s. the elongation of the spring in terms of the generalized coordinates u. The spring component acts as an energy storage device.

Another important component in structural mechanics is the viscous-damping element. Unlike the spring, the damping element dissipates energy from a vibrating structure. The most basic analytical model of a damping element is given by the linear relationship

f_d= c ˙e, (43)

where c is the coefficient of viscous damping and ˙e = ˙u₂− ˙u₁ is the relative velocity.

Figure 6 shows an illustration of a simple system with springs k and dampers c connected between two masses m1and m2. External forces F are applied to the bodies and u1and u2are the displacements of the respective body. Using Newton’s second Law it can be shown that the equation of motion for the system is

m1 0 0 m2

¨u1

¨ u2

+c1+ c2 −c2

−c2 c2

˙u1

˙ u2

+k1+ k2 −k2

−k2 k2

u1

u2

=F1(t) F2(t)

, (44)

which in matrix notation can be written as

M¨u + C ˙u + Ku = F(t). (45)

This equation is the general form of the equations of motion for a system with n-DoF and the equation can be claimed to be the fundamental equation in structural dynamics and linear vibration theory [8].

The coefficients M, C and K are called mass matrix, damping matrix and stiffness matrix, respectively, and all have the dimensions n × n. The vector u is called general displacements vector and is either expressed in physical or generalized coordinates and F(t) is the load vector, both have the dimensions n × 1. The off-diagonal terms in the stiffness matrix represents stiffness coupling.

The strain energy of the system, i.e. the potential energy due to the elasticity of the system, is given by [8, p. 228-238, 440–441]

U = 1

2u^TKu, (46)

and the kinetic energy of the system is given by T = 1

2u˙^TM ˙u. (47)

The damping matrix is in practice very difficult to determine. A significant problem is that the damping properties of a structure are frequency dependent. [4] Therefore, the damping is often formulated in terms of the mass and stiffness matrix, which is called Raleigh Damping. Another frequently used method when the system is described in modal terms is to assume the damping matrix to be diagonal and insert the damping as proportional factors on the diagonal elements. This is commonly referred to as modal damping. [8, p. 302-304] In the scope of this thesis modal damping with a damping factor of 3 % will be used in modeling of flexible bodies in SID format (see section 4.5). In section 4.6 a spring-damper model will be described, and for that model the damping parameters will be selected arbitrarily.

(18)

Figure 6: A 2-DoF spring-damper system.

3.4.2 Finite element method

The finite element method (FEM) is an efficient and commonly used method in structural dynamics. The method is usually implemented in finite element software. The method is in short based on the description of a structure as an idealized assemblage of individual structural elements. The basic idea being that it is easier to solve the equations that model the individual structural elements, and assembling these to describe the full structure, than trying the find a solution for the complete structure directly. The structural elements can take on a wide variety of shapes. Figure 7 shows an illustration of a structure, which have been divided into quadratic 2D elements. The dots represents the nodal points at which the elements are interconnected. In finite element analysis the structure is completely covered by these elements, in 2D or 3D depending on the analysis, so as no gaps are left between element domains. The procedure of dividing a structure into a system of discrete structural elements is commonly referred to as mesh generation. [4]

Figure 7: Quadratic 2D mesh. The black dots represent nodal points, which are associated with a specific number of DoF.

Using the finite element approach it possible to create a description of the full structure from the individual structural elements. The governing equation in finite element method for linear analysis in structural mechanics is equation 45. [4] The dimensions of the system, i.e. the number of DoF, is defined by the number of nodal points and the DoF of the respective nodal points. A full finite element model with a very fine mesh can therefore be represented by a system of thousands and thousands of DoF, resulting in very large system matrices (M, K and C). This results in systems that often demand very high computational power. An important subject is therefore model order reduction, which can be used to find representations of the system with less DoF. This will be further developed in section 3.5.

3.4.3 Beam theory

Considering the solid and homogeneous cantilever beam presented in figure 8. One side of the beam is clamped and the other end is completely free. A force and torque, respectively, is applied to the beam at a certain distance from the clamped end. The analytical expressions describing the deflection and rotation for a force in y direction are [33]

(19)

θ_x(α) = F_yl² 2EIx

β², δ_y(α) = F_yl³ 3EIx

β³ (48)

and for a torque in x direction

θx(α) = M_xl EIx

β, δy(α) = M_xl² 2EIx

β², (49)

where α + β = 1, l is the length of the beam, E is the Young’s modulus and Ixis the second moment of area.

Similarly, the deflection and rotation for a force in x direction is given by θ_y(α) = F_xl²

2EIy

β², δ_x(α) = F_xl³ 3EIy

β³ (50)

and for a torque around y

θy(α) = Myl EIy

β, δx(α) = Myl² 2EIy

β². (51)

The torsion (twist) of a beam can, by using Saint-Venants torsion theory, be expressed as [33]

θz=Mzl

GK, (52)

where K for a rectangular cross-section b × h (b > h) is given by K =bh³

3

F (b/h).¯ (53)

The factor ¯F (b/h) is a function, the value of which can be found tabulated in literature (see for example [33, p. 76] or [20, p. 334]) and G is the Shear modulus. There exists a relationship between Young’s modulus, the Shear modulus and the so called Poisson’s ratio v [20, p. 205]

G = E

2(1 + v). (54)

The elongation (stretch) of a beam can be calculated with δz= Fzl

EA, (55)

where A is the cross-sectional area of the beam.

The second moment of area is defined [20, p. 78]

I_x= Z

A

y²dA = bh³

12, (56)

which similarly for y is

I_y= Z

A

x²dA = b³h

12. (57)

Figure 8: Bending beam.

(20)

3.5 Model order reduction

Simulations of flexible multibody models are often dependent on the use of finite element models, where a large number of DoF are used to describe the deformation of a body. To enable simulations of large and complex systems it is therefore often necessary to reduce the size of the FE models. This is called model order reduction and will be presented in this section. A brief introduction of the general concept will be given directly below. After that, important subjects within the area will be presented in more detail.

As presented in section 3.4.2 the general form of the equations of motion of a structure is

M¨u + C ˙u + Ku = f , (58)

where the coefficient matrices M, C, K are the mass matrix, viscous damping matrix and stiffness matrix respectively, all with the size n × n. The displacement vector u, represents either physical or generalized displacements, and f is the load vector, and both vectors have the dimensions n × 1. A dot denotes differentiation with respect to time, t. In general, the objective of model order reduction is to find a reduced model with m DoF, where m << n, which still preserves essential characteristics of the full model. The general approach to reduce the order of a system is to apply a transformation on the displacement vector [10, 8]:

u = TuR, (59)

where T is a transformation matrix with dim(T) = n × m and uR is a reduced displacement vector with dim(uR) = m × 1. If the transformation is applied on equation 58 we obtain a reduced order system of the equations of motion

MR¨uR+ CRu˙R+ KRuR= fR, (60) where

M_R= T^TMT, C_R= T^TCT, K_R= T^TKT, f_R= T^Tf . (61) The coefficients MR, CR, KRare the reduced mass, damping and stiffness matrices, respectively, with the dimensions m × m and fR is the reduced load vector with dim(fR) = m × 1.

3.5.1 Floating frame of reference

Several methods on how to include both the large nonlinear motion and the deformations of a body in the modeling of multibody systems have been proposed in literature. Two efficient and commonly used methods are the absolute nodal coordinate formulation and the floating frame formulation [25].

The absolute nodal coordinate formulation [29, p. 309ff.] can be used in the analysis of flexible bodies that undergoes large deformations. For systems where the deformation of the body is small compared to the rigid body motion, the floating frame of reference [29, p. 189ff.] is preferably used. The floating frame of reference formulation is one of the most commonly used approaches of representing flexible bodies in computer simulations of multibody systems. The concept of this formulation is basically to describe the flexible body by using two sets of coordinates. One set of coordinates describe the location and orientation of a body reference coordinate system and another set of elastic coordinates describes the deformation of the body in terms of the body coordinate system. The origin of the body frame of reference do not need to be rigidly attached to the flexible body. It is however required that the there is no rigid body motion between the flexible body and the body frame of reference. It should be observed that the floating frame of reference formulation do not result in a separation of the rigid body motion and the elastic deformation [29, p. 191.]. Several different coordinate systems can be selected as the floating frame of reference, meaning that the motion of the body reference system do not necessarily coincide with the rigid body motion.

The floating frame of reference formulation is presented in figure 9, where Pi is a specific material point on a single flexible body i. The floating frame of reference formulation separates the description of the coordinates of point Piinto the coordinates of the body frame of reference Oi, i.e. the floating frame of reference, superimposed with the coordinates of the elastic deformation ¯uif expressed in reference frame O_i. This results in the following description of point P_i

rP = Ri+ Aiu¯i= Ri+ Ai(¯uiP+ ¯uif), (62)

(21)

where ¯uiP is the position of Piin the undeformed state and ¯uif is the deformation vector describing the displacement of Pi due to the deformation of the body. Both vectors are expressed in the floating frame of reference O_i. If the body is undeformed, the deformation vector vanish and the motion described by equation 62 results in the exact modeling of rigid body motion, i.e. equation 1. Using the floating frame of reference approach, it is possible to derive expressions for the velocity, acceleration and the equations motions for a flexible body. An extensive derivation of these expressions can be found in [29, p. 189ff.].

Figure 9: Floating frame of reference.

A strength of using the floating frame of reference approach in simulations of flexible multibody systems is the possibility of splitting the modeling of the system into two distinct parts. One side of the floating frame of reference approach is concerned with the modeling of the multibody dynamics, with important aspects as how to describe the kinematics and principles of mechanics of the multibody system. The other side is concerned with the description of the deformation of the different flexible parts of the system, where starting from a finite element formulation of a body, the FE model can be put through a model order reduction process to enable faster simulations. [25]

The following subsections are concerned with the concepts dynamic substructuring and model order reduction.

3.5.2 Dynamic substructuring

The concept of dynamic substructuring is concerned with ways to capture the dynamics of large and possibly complex structures by dividing the full structure into a number of smaller substructures or components, under the assumptions that it is easier to determine the dynamic behavior of each substructure respectively than for the complete structure. The dynamics of the complete structure is then evaluated by assembling the substructures. [34] Dynamic substructuring has been found to be very efficient when the full structure consists of several natural substructures or components. [8] Applications of this kind of natural substructuring can be found in the automobile industry, where components such as the body and frame constitutes a natural division of the larger system. A study where dynamic substructuring is applied on a passenger car can be found in [35]. Another structure where it is possible to find several natural substructures are industrial robots. In figure 1 these substructures constitute the different parts:

base, frame, lower arm, armhousing and so on.

Following the introduction given above, the concept of dynamic substructuring is intuitively clear:

divide a complex structure into substructures, which can be assumed to be easier to analyse. To efficiently accomplish this, it is however necessary to take some considerations into account. Therefore, a brief overview of the classification of different substructuring approaches and some theoretical considerations will be presented below [17, 34]. Dynamic substructuring can broadly be divided into three domains:

• The physical domain, in which the structural properties of a substructure are described in terms

(22)

of the mass, stiffness and damping properties of the substructure. These are given by the stiffness, mass and damping matrices of the substructure, i.e. the full FE model.

• The modal domain, where the dynamic behavior of the substructure is described by its modal responses.

• The frequency domain, where the substructure is interpreted in terms of its frequency response functions.

It should be observed that the three domains noted above, at least from a theoretical perspective, for a specific substructure will contain the same information. This under the assumption that no model order reduction process is applied on the structure. [17]

Since dynamic substructuring is based on the concept of dismantling a full structure to smaller substructures, only to later reassemble the full structure, an important question is the coupling of two or more substructures. Regardless of the specific method used for the coupling, two fundamental conditions need to be satisfied [17]:

1. Compatibility condition: the requirement that the substructures need to be compatible in regards to the displacements at the interfaces, i.e. the DoF of the displacements of the connecting interfaces of the substructures need to be same.

2. Equilibrium condition: the forces of the coupled substructures’ interface DoF must be in equilibrium, i.e. the forces need to be opposite in direction and equal in magnitude.

If the dynamic models of the substructures are known, the response of the coupled system can be calculated if the coupling conditions are satisfied. In dynamic substructuring it is possible to choose displacements or forces as unknowns at the interface. It should, however, be noted that some coupling methods choose both interface displacements and forces as unknowns.

3.5.3 Component-mode synthesis

When studying large and complex structural dynamics problems, it is often necessary to use some form of component-mode synthesis (CMS) method [8, p. 532ff.]. CMS methods usually refers to dynamic substructuring where the substructure is a reduced space and the term ”mode” should be understood as any vector representing the reduction space, e.g. a Ritz vector, which is a basis vector used to describe the displacements of points within a substructure. Usually, CMS methods includes a modal representation of the substructure, from which modal vectors are obtained. These modal vectors are then used to reduce the equations of motion of the physical system, reducing the system to the modal domain. By performing this reduction, the full set of physical coordinates are reduced to a smaller set of generalized coordinates.

[17]

The equation of motion of an undamped elastic substructure can be expressed as

M^su¨^s+ K^su^s= f^s+ r^s, (63)

where the superscript s designates the specific substructure, M^sand K^sare the substructure’s mass and stiffness matrix, respectively, with the dimensions n × n. The displacement vector u^shas the dimensions n × 1. The coefficient matrices and the displacement vector are expressed in the substructures original physical coordinates. The forces on the right hand side of equation 63 consists of the external forces applied to the substructure f^s and the reaction forces r^s on the substructure due to the coupling to adjacent substructures at the interface DoF. The damping of the system is in CMS methods usually treated as modal damping imposed on the modes of the reduced-order system, and is therefore omitted from equation 63 (see section 3.4.1). Some special CMS methods have been proposed with viscous damping included in the model but these are not in wide use. [8, p. 532ff.] For the purpose of this thesis the damping will be imposed on the modes of the reduced-order system.

The main idea in CMS is to express the physical coordinates u in terms of the component generalized coordinates ˜u with dim(˜u) = m × 1, by using a coordinate transformation

u^s= Ψ^su˜^s, (64)

(23)

where the transformation matrix (often also called component-mode matrix ) Ψ^s with the dimensions n × m represents the reduction basis and contains either the preselected component (assumed) modes or the Ritz-vectors. The main difference between different CMS methods is usually the choice of these assumed modes (see section 3.5.4 for a brief overview of common modes used in CMS methods). For the reduction to be efficient it is usually required that the reduction basis significantly reduces the DoF of the substructure (m << n). If the coordinate transformation in equation 64 is used in equation 63 the equation of motion of the substructure s is reduced to

M˜^su¨˜^s+ ˜K^su˜^s= ˜f^s+ ˜r^s, (65) where

M˜^s= Ψ^sTM^sΨ^s, (66)

K˜^s= Ψ^sTKΨ^s, (67)

˜f^s= Ψ^sTf^s, (68)

˜r^s= Ψ^sTr^s. (69)

The coefficients ˜M^s and ˜K^s, with the reduced sizes m × m, are the reduced mass and stiffness matrices, respectively, and ˜f^s and ˜r^sare the reduced force vectors, both with the dimensions m × 1.

Several CMS methods have been proposed in literature [8, 17]. Two of the most common methods are Guyan reduction and the Craig-Bampton method. These two methods are presented in section 3.5.5 and 3.5.6, respectively.

3.5.4 Modes used in CMS

As observed in section 3.5.3 the main difference between different CMS methods are the choice of assumed modes as the reduction basis. All kinds of modes can be selected as the reduction basis, e.g. eigenmodes, normal modes, attachment modes etc. This section will give a brief overview of some of the most common modes used in CMS methods.

To derive different types of component modes to be used as a reduction basis in the transformation matrix Ψ^s in equation 64 it is useful to divide the DoF of the structure into internal DoF, denoted by the subscript i, and boundary DoF, denoted by the subscript b [8, p. 534ff.]. The internal DoF are the DoF that will be condensed and the boundary DoF are the remaining DoF that will be retained in the reduction process. It is also common practice to refer to the reduced internal DoF as the slave DoF and the retained boundary DoF as the master DoF. The equation of motion in equation 63 can be partitioned into boundary DoF and internal DoF as (the superscript s will be omitted to increase readability)

M_ii Mib

M_bi M_bb

¨ui

¨ u_b

+K_ii Kib

K_bi K_bb

u_i u_b

=

0i

f_b+ r_b

, (70)

where the boundary DoF are kept in ub and the internal DoF are kept in ui. The internal force fi is omitted (assumed to be zero), since it assumed that if an internal DoF is exposed to a load it will be relabel as a boundary DoF [8, p. 533-534].

Normal modes are eigenvectors and can be classified according to the interface boundary conditions for the specific substructure: fixed-interface normal modes, free interface normal modes, hybrid-interface normal modes, or loaded-interface normal modes.[8, p. 534-535].

Fixed-interface normal modes can be obtained by constraining all the retained DoF and evaluat- ing the eigenmodes of the substructure by solving an eigenvalue problem. Starting from the partitioned system in equation 70 all the retained DoF are constrained, i.e. ub = 0, which results in

Mii¨ui+ Kiiui= 0. (71)

The system in equation 71 can be solved with the following eigenvalue problem [8, p. 535]:

(Kii− ω_i,j² Mii)φ_i,j= 0, j = 1, 2..., Ni (72)

Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

Examensarbete 30 hp Juni 2019

Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

Christoffer Fjellstedt

Abstract

Methods for including stiffness parameters from reduced finite element models in simulations of multibody systems

Christoffer Fjellstedt

1 Popul¨ arvetenskaplig sammanfattning

Contents

2 Introduction

2.1 Goals of the thesis

2.2 About the project

3 Theory

3.1 Modelling of robot manipulators

3.2 Multibody systems

3.3 Constrained multibodies

3.4 Structural mechanics

3.5 Model order reduction