Safe trajectory of a piece moved by a robot

ESGI 158 Report Safe trajectory of a piece moved by a robot November 9, 2020 Ernest Benedito, ernest.b[email protected], Oliver Bond, [email protected], Thomas Babb, [email protected]x.ac.uk, Juan R. Pacha, [email protected], Sandeep Kumar, [email protected], Joan Sol`a-Morales, [email protected]. 1 Introduction The company F.EE (www.fee.de/en.html) is an international company supplying automation technologies for industry, mostly in the automotive sector. The company is interested in determining the equations of the trajectories and the orientations of pieces that are moved by the action of a robot with seven degrees of freedom, of which six correspond to rotations of the arms and one to longitudinal transfers. In many cases, the piece is a curved metal sheet, and in the others, a three-dimensional object. At the initial moment, the piece is at a certain point/location, orientation, and rest position, then the robot grabs and leaves it to another given point/location, orientation, and rest position corresponding to the other time instant. During its movement, the piece can suffer reversible and irreversible deformations caused by both mass forces and surface forces. The mass forces are due to translation and rotations (depending on accelerations and angular velocity of the body) while surface forces are due to air drag (depending on the velocity of the body) that act on it. The company wants that the resulting motion does not result in irreversible deformations on the piece and that the time between the initial and final moments does not exceed a certain threshold. The actual movement of a piece is a consequence of the movements of the robot. This in turn depends on the design of the trajectory to be followed by the piece, as well as on the robot’s ability to faithfully follow the desired trajectory. Any solution to the problem must consider not only the conditions under which the desired trajectory cannot cause irreversible deformation but also, the conditions under which the robot can follow the trajectory. The former depends on the desired velocity, acceleration, angular velocity of the piece, and the latter depends on the so-called jerk, or, the first-derivative of the accelerations. Certainly, the problem is complex and is rarely addressed in the scientific literature. There is abundant literature on the trajectory design of a piece moved by a robot with varying degrees of freedom [1, 2] but the deformations caused by mass forces or friction are generally not taken into account. This is certainly a matter of practical interest to the manufacturing industry and that it should be for the scientific community as well. 1 arXiv:2011.03330v1 [cs.RO] 22 Oct 2020 Figure 1: Left: Metal piece moved by the arm of the robot. Right: Metal piece as a part of a car. Due to its complexity, we have understood that the problem cannot be solved in limited time such as five days of the ESGI meeting, but in this time, a conceptual mathematical framework has been laid out, which can be later adapted and solved. Therefore, during the ESGI days, we have restricted the objectives of the problem to the following tasks: •Perform a literature search in this area, including a search of numerical codes to solve these kinds of problems. •Write the equations of motion of the body-fluid system, in the most general case possible, by determining the stresses and deformations suffered by a body due to the forces of mass and friction acting on it when being moved by another. •Together with the possible use of commercial, discuss numerical algorithms to solve these equations. •Define a very simple case of the problem. This has been the movement of a longitudinal piece, held at one end, to a body that can move in a direction perpendicular to that of the force of gravity, and can rotate the piece in the plane defined by the direction of motion of the body and that of the force of gravity. The results of the first point will be described in Section 5 below, together with conclusions about possible future work, and the next three points will be described in Sections 2, 3, and 4. We want to point out that we have devoted our efforts to the so-called direct problem, rather that the more difficult, but perhaps more important in practice, inverse problem. Let us explain the difference: In the direct problem, one assumes that the trajectory of the arm and its rotations are known and given, and the problem is to calculate deformations and stresses of the piece. In the inverse problem, the goal is to determine trajectories and rotations such that the induced deformations and stresses satisfy the restriction requirements given by the nature of the pieces, to avoid permanent deformations. Even if one considers that only the inverse problem is of practical interest, it is clear that the solution of the inverse problem requires the solution of many direct problems, to achieve the optimal candidate. 2 The equations Before explaining the equations and boundary conditions specific to plate theory, we provide a brief introduction to the governing equations in solid mechanics. This explanation is a condensed version of that found in [3], which gives a much more rigorous introduction to the field. 2 Introduction to solid mechanics Suppose that any point inside an elastic object P ⊂ R3has an initial Lagrangian coordinate ~ X∈ P, which is fixed in the material. At time t > 0, this coordinate changes with elastic deformations of the object, into an Eulerian coordinate ~x(~ X, t), where the coordinate ~x remains fixed in space. From these, we can construct the displacement field ~u(~x, t) := ~x(~ X, t)−~ X, which in part is found by formulating an equation relating body and surface forces, and using conservation of momentum to obtain a partial differential equation known as the Navier equation in which ~u(~x, t) is the dependent variable. Consider an arbitrary volume Vcontained within the elastic material. This small volume experiences body forces (e.g. due to gravity) and surface forces. In the case where the density ρis uniform throughout the material, one can write down expressions for both of these effects, as well as the rate of change of momentum (where momentum is the mass of an object multiplied by its velocity). In this case, the ith component of the displacement field, ui, satisfies d dtZZZV ∂ui ∂t ρdP | {z } rate of change of momentum =ZZZV giρdP | {z } total body force +ZZ∂V σijnjdS | {z } total surface force , where giis the i-th component of the body force, and σij is the (symmetric) stress tensor, which tells us the ith component of the force per unit area acting at a point on the surface with outward normal nj. By taking the time derivative inside the first term (using the fact that the density is constant) and applying the divergence theorem on the final term, we can use the fact that Vis arbitrary and the assumption that each integrand is continuous to obtain Cauchy’s momentum equation ρ∂2ui ∂t2=ρgi+∂σij ∂xj .(1) To be able to determine the displacement field relies on us knowing the stress tensor σij. Thankfully, there is a constitutive relation, known as Hooke’s law, which postulates a linear relationship between the stress and another quantity called the strain, a dimensionless quantity signifying the elastic extension of an object relative to its original state. To identify the exact form of the strain tensor, we can consider two particles with positions ~ Xand ~ X+δ~ X, displaced to ~x =~ X+~u(~ X, t) and ~x +δ~x =~ X+δ~ X+~u(~ X+δ~ X, t). Taylor’s theorem can be used to show that (upon neglecting quadratic terms) |δ~x|2=δ~ X+ (δ~ X· ∇ ~ X)~u(~ X, t) 2, and therefore that |δ~x|2−δ~ X 2= 2εijδXiδXj, where the strain tensor εij is given by εij =1 2∂ui ∂Xj +∂uj ∂Xi +∂uk ∂Xi ∂uk ∂Xj∼1 2∂ui ∂Xj +∂uj ∂Xi.(2) If the stress and strain are scalars, then Hooke’s law simply states that σ=Eε, where Eis a constant known as the Young modulus. However, in this context, σij and εij are both rank-2 tensors, so a linear relation between them must involve a rank-4 tensor, so that σij =Cijklεkl, i, j, k, l = 1,2,3. Here, Cijkl has 81 entries in total. However, by considering the stress acting at the surface of the volume V(by considering a pillbox-type argument), one can show that the rank-2 tensors are symmetric (i.e. 3 Figure 2: A simple illustration of the coordinate system used in Kirchoff-Love plate theory for a (flat) plate. σij =σji and εij =εji). Under further modelling assumptions, such as homogeneity and isotropy of the material, one can derive a stress-strain relation of the form σij =λεkkδij + 2µεij,(3) where λis the bulk modulus and µis the shear modulus of the material. Both of these constants are collectively called the Lam´e constants and they are material properties which indicate a material’s tendency to withstand deformation. Upon substituting (3) and (2) into (1), one obtains the Navier equation, written in vector form as ρ∂2~u ∂t2=ρ~g + (λ+µ)∇(∇ · ~u) + µ∇2~u. (4) By imposing suitable boundary conditions, one can solve equation (4) for the displacement field ~u(~ X, t). The usual techniques of applied mathematics, such as asymptotic analysis and numerical methods, can be used to solve this equation, although numerical methods will be much more appropriate for arbitrary geometries (such as those from CAD files provided by component manufacturers). Classical plate theory A considerable simplification in our situation is that many components moved around by a robot are three-dimensional but are thin; for example, a car door. Although many people think of this type of component as strong, namely due to being made of metal, thin metal sheets are prone to elastic displacements which could exceed an elastic limit, deforming plastically and then being unsuitable for use in any further manufacturing (and thus needing to be recycled). Thankfully, a theory has been developed to tackle the case where the component is large and thin. This is called Kirschoff-Love theory (or classical plate theory), where in-plane displacements are disregarded, and that displacements normal to the plane are considered relative to a mid-plane: a surface equidistant between the top and bottom face of the plane. This idea is illustrated in Figure 4, although we emphasise that the plate need not be uniform everywhere since the coordinates are local to the plate. We guide the reader through the basics of this theory; and refer them to [4] for a much more detailed explanation. In the specific case shown in Figure 4, an educated guess for the displacement field is given. This is of the form ui(x, y, z, t) = X j (z)jφ(j) i(x, y, t), 4 where φ(j)are functions chosen so that the principle of virtual displacements is satisfied and zis the coordinate in the direction of thickness. A virtual displacement of an elastic system is an infinitesimal displacement which can one might suspect to arise based on the current configuration of forces and if it was perturbed slightly. For example, if an elastic beam of length Lis fixed at a wall at x= 0 and being pulled away by a force Fbeing applied at its free end, then the beam can be thought of as having a “virtual” displacement of δu(L) at its free end arising from a “virtual” force δF. The idea of the displacements being “virtual” is so-called because they are fictional and are unrelated to the displacements and actual loads on the system. Virtual work and virtual displacement The idea of “virtual work” arises in two possible ways. We firstly recall that the work done by a force is the product of its projection (dot product) onto its direction of displacement, and its magnitude (or, put simply, “work is force times distance”). Virtual work can be considered to arise either from (a) an actual force moving through a virtual displacement, or (b) a virtual force moving through an actual displacement. If a lowercase delta before a variable denotes an infinitesimal virtual change in that variable, then the virtual work is given by δW =ZZZP ~ F·δ~udP. Virtual work can be classified as either external or internal. Forces applied by external sources will do work when moving through virtual displacements, and to this effect one can write down the total external virtual work. When an elastic body is subjected to body forces of ~ fper unit volume and surface tractions ~ Tper unit surface area (where Γσis the subset of ∂Pon which stresses are specified), upon moving through virtual displacements δ~u it will do a total amount of virtual work δV due to applied forces. This is given by δV =−ZZZP ~ f·δ~udP+ZZΓσ ~ T·δ~udS.(5) A deforming elastic object will also undergo internal stresses and then the work can be obtained in terms of the total strain inside it rather than due to specific vector forces. In simple terms, the strain of an elastic body is the ratio of the extension of an elastic body to its original size, and is therefore a dimensionless quantity, and since it involves derivatives of displacement, virtual strains δεij can be expressed in terms of virtual displacements δu. By considering work done by normal stresses and shear stresses through virtual displacements, one can show that the total internal work δU is given by δU = 3 X i=1 3 X j=1 ZZZP σijδεijdPwhere δεij =1 2∂δui ∂xj +∂δuj ∂xi.(6) At this point, we can concisely state the principle of virtual displacements, i.e. that a continuous body in equilibrium will have a total virtual work by actual forces through virtual displacements of zero. Put simply, δV +δU = 0; as each of these terms are integrals, the problem then becomes about minimising these integrals (through a weak formulation) rather than solving for an unknown function with a differential equation (through a strong formulation). More specifically of interest to us, the dynamical version of the principle of virtual displacements is ZT 0 (δU +δV −δK) dt= 0,(7) where δK is the virtual kinetic energy of the system. This method of deriving governing equations for an elastic body is an attractive alternative to equation (4) because (a) it does not need to rely on the use of constitutive laws upfront (e.g. Hooke’s law), and (b) the elasticity problem can be posed 5 in the form of a variational problem, thereby making it amenable to both analytical treatments (e.g. calculus of variations, Euler-Lagrange equations, Hamilton’s principle) and contemporary numerical methods such as finite element methods (FEM). How one might approach this problem using FEM is to be discussed later. Deriving the governing equation We are interested in Kirchoff’s hypothesis which provides the elastic displacement field a priori as u(x, y, z, t) = u0(x, y, t)−z∂w0 ∂x , v(x, y, z, t) = v0(x, y, t)−z∂w0 ∂y , w(x, y, z, t) = w0(x, y, t), where u0,v0and w0are displacements from the mid-plane in the x,yand z-directions respectively. In particular, we can set u0≡0 and v0≡0, due to neglecting in-plane displacements. Namely, we are interested in the domain [−w, w]× F, where Frepresents the cross-section of the plate at z= 0. Substituting the above form of the displacement field into the linearised strains (2), neglecting in-plane displacements, it is straightforward to show that the strain tensor is given by 11 =1 2∂w0 ∂x 2 −z∂2w0 ∂x2 22 =1 2∂w0 ∂y 2 −z∂2w0 ∂y2 12 =1 2∂w0 ∂x ∂w0 ∂y −2z∂2w0 ∂x∂y  13 =1 2−∂w0 ∂x +∂w0 ∂x = 0 23 =1 2−∂w0 ∂y +∂w0 ∂y = 0 33 = 0, (8) where we have given six entries instead of nine, owing to the fact that the strain tensor is symmetric (so that ij =ji). This puts us in a position to calculate the internal virtual work δU, given by δU =Zw −wZF (σ11δ11 + 2σ12δ12 +σ22δ22) dFdz =−Zw −wZFσ11z∂2δw0 ∂x2+ 2σ12z∂2δw0 ∂x∂y +σ22z∂2δw0 ∂y2dFdz+ h.o.t. =−ZFM11 ∂2δw0 ∂x2+ 2M12 ∂2δw0 ∂x∂y +M22 ∂2δw0 ∂y2dF+ h.o.t. where “h.o.t.” is an abbreviation for “higher-order terms”; terms which are quadratic in the partial derivatives of δw0and can therefore be neglected (although this is no longer appropriate if von Karman strains are considered). We also have Mij =Zw −w σijzdz, (9) are the stress moment resultants. We emphasise that in general, there will also be stress resultants Nij =Zw −w σijdz, 6 but neglecting the terms which are quadratic in the derivatives of w0leads to them being absent from the expression for the internal virtual work. Meanwhile, the external virtual work is given by δV =−Zw −wZF ~ f·δ~udFdz−Zw −wZ∂F ~ T·δ~udSdz(10) =−Zw −wZF (q−kw0)δw0dxdy. (11) The second term in equation (10) disappears because the traction force, Ti=σijnjwhere ~n is the normal to ∂F(a closed curve in the (x, y)-plane), is perpendicular to the z-direction in which the displacement is solely assumed to take place. In equation (11), q(x, y) is the net load on F, and the −kw0contribution arises from Hooke’s law, where kis the stiffness constant of the plate material. The internal kinetic energy δK is given by δK =ZPZw −w ρ( ˙uδ ˙u+ ˙wδ ˙w+ ˙wδ ˙w) dzdxdy =−I0ZP ( ˙u0δ˙u0+ ˙w0δ˙w0+ ˙w0δ˙w0) dxdy−I2ZP∂˙w0 ∂x ∂δ ˙w0 ∂x +∂˙w0 ∂y ∂δ ˙w0 ∂y dxdy(12) where the moments of inertia are given by Ik=Zw −w zkρdzwhere k= 0,1,2. We do not present the full derivation of the governing equations here due to the amount of algebra involved, but briefly describe how it is done (see pages 103-105 of [4] for a full derivation). By substituting results (6), (11) and (12) into (7), making use of the virtual strains, employing techniques from the calculus of variations and setting the coefficient of δw0to zero, one can show that ∂2M11 ∂x2+ 2∂2M12 ∂x∂y +∂2M22 ∂y2−kw0+q=I0 ∂2w0 ∂t2−I2 ∂2 ∂t2∂2w0 ∂x2+∂2w0 ∂y2. Finally, it is desirable to express this governing equation in terms of displacements rather than moments. For a homogeneous, isotropic plate, a constitutive law tells us that   σ11 σ22 σ12  =E 1−ν2  1ν0 ν1 0 0 0 1 −ν    11 22 12  ,(13) where Eis the Young modulus of the material, and νis the Poisson ratio. Combining the constitutive law (13) with (9) and the strain-displacement relations (8), one obtains D∇2∇2w0=−q(x, y, t)−2ρh∂2w0 ∂t2,(14) where the bending stiffness,D, is given by D=2h3E 3(1 −ν2),(15) which will be loosely referred to as the “Kirchoff-Love equation”. 7 Initial and boundary conditions Along with the equation of motion (14), several boundary conditions based on likely physics in a factory setting need to be imposed. Let C ⊂ ∂Fbe the part of the large face boundary which is clamped by the robot, then the plate 1. is initially undeformed, i.e., w0(x, y, t = 0) = 0, 2. is initially stationary, i.e., ∂w0 ∂t (x, y, t = 0) = 0, 3. does not deform where it is clamped, i.e., w0(x, y, t) = 0, x, y ∈ C; 4. and no bending moments or loads at the free boundary: ∇2w0(x, y, t)=0,∂ ∂n(∇2w0(x, y, t)) = 0,(x, y)∈∂F\C Changing the frame As it stands, it suffices to solve equation (14) for the transverse elastic displacements of the plate, provided that the plate is not being moved externally. Considering how the plate is being moved around by a robot, this is not satisfactory; if the plate is being moved in a non-inertial frame, it will experience “fictitious” forces, namely the Euler, Coriolis and centrifugal forces. These forces are so-called because they arise from a change of frame rather than a physical mechanism. They are mathematical in nature and should therefore be incorporated within the external load q(x, y, t) as it appears in equation (14). The mathematical techniques for changing frames are standard and can be found in classical mechanics texts such as [5]. We use ˆ Sto denote an inertial frame which remains fixed over time, with origin ˆ Owhich is fixed in space, and orthonormal basis {ˆ ~e1,ˆ ~e2,ˆ ~e3}. We use Sto denote the non-inertial frame, which is fixed local to the plate and whose origin Omoves with angular velocity ~ω relative to ˆ O. The relationship between the two coordinate systems is illustrated in Figure 3. F.EE have complete control over the rotational movement of the robotic arm, as well as translational motion along a straight line. In practice, the rotational motion of the robot motion may be prescribed by a rotation matrix Rij(t) which acts as a time-dependent linear transformation from ˆ S to Sand may be expressed in terms of Euler angles. If ˆ D and D denote time derivatives in the inertial and non-inertial frames respectively, and ~r is the position vector of a point in Sfixed on the plate relative to O, then the Coriolis formula is given by ˆ D~r = D~r +~ω ×~r, where ~ω is the angular velocity vector whose entries satisfy Rij ˙ Rji = 3 X k=1 ijkωk.(16) Applying the Coriolis formula twice can be used to obtain an expression for the acceleration ˆ ~a in the inertial frame in terms of the acceleration ~a in the non-inertial frames: ˆ ~a =~a + (D~ω)×~r + 2~ω ×D~r +~ω ×(~ω ×~r) + ~ A, (17) where ~ A=ˆ D2~x is the acceleration of Orelative to ˆ S. 8 Figure 3: An illustration of plate motion in two frames of reference. This is useful because it can be used in principle to re-express the Kirchoff-Love equation (14) in a non-inertial frame local to the plate, as the acceleration arises naturally there. The acceleration formula (17) can then be used to to re-express Newton’s second law in Srather than ˆ S: mˆ ~a =~ Fin ˆ S=⇒m~a =~ F−m(D~ω)×~r | {z } “Euler force” −2m~ω ×D~r | {z } “Coriolis force” −m~ω ×(~ω ×~r) | {z } “Centrifugal force” +m~ Ain S.(18) Since we are only concerned with elastic oscillations of the plate in the ~e3-direction, it suffices to take the dot product of equation (18) with ~e3in order to obtain the equivalent of the Kirschoff-Love equation (14). We do not write down the equation in full, since it depends on the specific rotations being applied. However, we simply state that the term involving the second derivative of w0with respect to time corresponds to a·e3and that the external loads q(x, y, t) correspond to ~ F. It is not quite enough to simply apply a change of frame to the Kirschoff-Love equation to take all physical considerations into account. There is also the weight of the plate, the inertia, the internal elastic forces inside the plate, and also the air drag. However, in a factory setting, the effects of these terms would be questionable, and may make the equations significantly more coupled without adding much insight. This is particularly the case for the air drag because a complete description would rely on coupling the elastic displacement equations discussed hitherto with the Navier-Stokes equations of fluid dynamics. A simplified starting point could be to introduce a simple drag law, where the drag force on an object in a fluid is proportional to the square of the speed at which object is being passed by the fluid. The exact nature of the physical effects is well beyond the scope of this report, but would be interesting nonetheless. Another physical phenomenon to consider is that many elastic materials deform plastically when their displacements pass a limit that is large enough in size. This is often called yield, and a hypothesis governing an instance in which a material deforms is known as a yield criterion. Yield hypotheses are frequently given in terms of functions of the entries of the stress tensor σ; for example, the von Mises’ 9 grant MTM PGC2018-100928-B-I00. S.K. is supported by the grant Severo Ochoa SEV-2017-0718. J.S.-M. acknowledges partial support by MINECO (Spain) grant MTM2017-84214-C2-1-P. References [1] D. E. Whitney, “The mathematics of coordinated control of prosthetic arms and manipulators,” 1972. [2] V. Braibant and M. Geradin, “Optimum path planning of robot arms,” Robotica, vol. 5, no. 4, pp. 323–331, 1987. [3] P. Howell, G. Kozyreff, and J. Ockendon, Applied solid mechanics. Cambridge University Press, 2009, vol. 43. [4] J. Reddy, Theory and Analysis of Elastic Plates and Shells, 2nd ed. CRC Press, 2007. [5] H. Goldstein, C. P. Poole, and J. L. Safko, “Classical mechanics third edition, chapter 7,” 2002. [6] R. v. Mises, “Mechanik der festen k¨orper im plastisch-deformablen zustand,” Nachrichten von der Gesellschaft der Wissenschaften zu G¨ottingen, Mathematisch-Physikalische Klasse, vol. 4, pp. 582–592, 1913. [7] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 5th ed. Oxford: Butterworth Heinemann, 2000, vol. 2. [8] A. M. Lush, “Coupling of beam and shell finite elements for the rapid analysis of tubular structures,” MsC Thesis, School of Engineering and Computing Sciences. Durham University, 2014. [9] K. Bathe, Finite Element Procedures. Prentice Hall, 2006. [10] A. Moysidis and V. Koumousis, “Hysteretic shell finite element,” Journal of Engineering Mechanics, vol. 145, no. 5, 2019. [11] Y. W. Kwon and H. Bang, The Finite Element Method using MATLAB, 2nd ed., ser. CRC Mechanical Engineering Series. Chapman & Hall/CRC, Boca Raton, FL, 2000. [12] Wikipedia contributors, “Ansys — Wikipedia, the free encyclopedia,” https://en.wikipedia.org/ w/index.php?title=Ansys&oldid=963235585, 2020, [Online; accessed 23-June-2020]. [13] ——, “ADINA — Wikipedia, the free encyclopedia,” https://en.wikipedia.org/w/index.php? title=ADINA&oldid=875541407, 2018, [Online; accessed 23-June-2020]. [14] ——, “FEniCS Project — Wikipedia, the free encyclopedia,” https://en.wikipedia.org/w/index. php?title=FEniCS Project&oldid=957113708, 2020, [Online; accessed 23-June-2020]. 16