An Interactive Physically-based Model for Active Suction Phenomenon Simulation

—While suction cups are widely used in Robotics, the literature is underdeveloped when it comes to the modelling and simulation of the suction phenomenon. In this paper, we present a novel physically-based approach to simulate the behavior of active suction cups. Our model relies on a novel formulation which assumes the pressure exerted on a suction cup during active control is based on constraint resolution. Our algorithmic implementation uses a classiﬁcation process to handle the contacts during the suction phenomenon of the suction cup on a surface. Then, we formulate a convenient way for coupling the pressure constraint with the multiple contact constraints. We propose an evaluation of our approach through a comparison with real data, showing the ability of our model to reproduce the behavior of suction cups. Our approach paves the way for improving the design as well as the control of robotic actuators based on suction cups such as vaccum grippers.


I. INTRODUCTION
Suction cups are widely used in Robotics to achieve interactions with nearby objects or to manipulate tools for interacting with the environment.Many different types of robots hold suction cups, such as traditional grippers, bioinspired robots [1], or soft robots [2], making them able to grasp an object or to climb along a surface.Robots equipped with vacuum grippers are used for various applications ranging from medicine to industry and manufacturing.Designing suction cups for robot manipulation purposes is crucial in many cases to ensure accurate control of their actions.However, few methods have been proposed as of today for simulating the behavior of suction cups (see Section II).Suction is a complex mechanical phenomenon, where both the physical soft properties of the suction cup as well as the properties of the surface on which the suction cup is stuck are of importance to perform a given task.The main bottleneck is the handling of deformable objects where the number of contacts between the surface and the suction cup is high, leading to a high computational cost.In this paper, we propose a novel physically-based model of active suction cup phenomenon using a constraint-based simulation (see Figure 1).Our deformable suction cup is simulated using the co-rotational Finite Element Method (FEM) [3].In our work, we focus on active phenomena where the pressure inside the suction cup is imposed by a regulator.The main contributions are the following: • a novel physics model based on a constraint formulation to simulate active suction cup behavior, and especially the pressure inside the suction cup cavity; • a classification algorithm of the nodes located on the inner surface of the suction cup for handling the contacts and the pressure exerted on the suction cup; • an evaluation of the model through a comparison of simulated data with real data.The remainder of the paper is organized as follows: Section II describes the related work on suction cup modeling, physically-based simulation of deformable objects and contact modeling.Then, Section III presents our model as well as its implementation.Section IV proposes the validation, both on simulated and real data while Section V ends the paper with conclusion and perspectives.

A. Suction Cup Modeling
There are relatively few studies about the modeling and simulation of the suction cup phenomenon in the litterature.One of the earlier but relevant works on the topic is the one from Bahr et al. [4] who designed and analyzed a climbing robot which was equipped with suction cups at the tip of its two legs.In their mathematical model, the authors formulated the contact phenomena using Lagrange constraints, but did not provide simulations of their model.Some years later, Liu et al. presented an analytical modelling of suction cups for window-cleaning robots [5], in which friction was omitted by considering that the window surface was always perfectly wet.Through a quantitative study of the attachment and detachment of a passive suction cup, Ge et al. modeled the forces acting on the cup considering only the vertical ones for the sake of simplicity [6].More Recently, Mahler et al. proposed a compliant suction contact model [7] wherein the suction cup was modeled as a quasi-static spring system.The model estimates the deformation energy to maintain a seal and quantifies the ability to resist external wrenches.Nevertheless, and to the best of our knowledge, neither physics-based model nor interactive simulations of suction cup phenomena exist in the literature.Such models require the use of interactive simulations of both a deformable model and a contact model for handling the contact between the deformable suction cup and different objects.

B. Physically-based Simulation of Deformable Objects and Multiple Contacts
The Finite Element Method (FEM) is widely used to obtain realistic physics-based simulations of deformable solids.Given the large displacements experienced by the suction cup during deformation, a purely linear FEM model cannot be used because it would lead to deformation artifacts.A corotationnal approach was proposed by Müller and Gross [8] to overcome this problem, in which the rotationnal part of the deformation is computed separately through a polar decomposition.Later improvements were also proposed by Muller et al. [3] allowing the interactive simulation of many deformable objects in a realistic way.
Concerning the simulation of suction cups, contact handling between deformable objects remains a bottleneck of the current simulations.It can be achieved with a constraint-based formulation, revealing a quadratic problem (QP) [9] [10] or a complementarity problem (CP).Quadratic methods build and solve a single system of constraints whereas complementary problem methods solve the contact problem in the constraint space: roughly, it means that each simulated object has its own system of constraints.The systems are firstly solved separetely, then the solver resolves the contact equations between adjacent objects.The first interactive contact model for deformable objects was proposed by Duriez et al. based on the Signorini's law [11].In [12], Otaduy et al. achieved contact handling by mixing a constrained dynamics formulation with semi-implicit nonpenetration constraints leading to a Mixed Linear Constraint Problem (MLCP).The high amount of contact points is often a bottleneck in real-time simulation as it leads to numerous contact constraints.Talvas et al. recently proposes an aggregate constraints formulation to tackle this issue [13].The method lumps contact points and consequently reduces the number of constraints.

A. Method Overview
Before describing our approach for suction cup simulation, we first formulate the dynamics of a physical system with generic constraints.According to Newton's second law, the dynamics of a body within the simulation follow: where M is the mass matrix, x the position of the degrees of freedom (i.e. the nodes of the FEM mesh), v = ẋ the vector of velocities.f provides the internal forces and f ext the external forces (due to contact and other constraints).Using implicit Euler integration with linearized forces, and given a time step h and current state (x 0 , v 0 ), the velocity increment dv on a time step is obtained by solving the linear system: We now consider our system as a deformable suction cup with a tetrahedral volumetric mesh and a triangular surface mesh whose vertices collide with triangles of another body.In order to prevent their interpenetrations, and to simulate contacts and pressure phenomenon, we solve the previous system for both bodies under contact and pressure constraints.These constraints provide the external forces f ext , which are defined as the product of H, a matrix of constraint directions, and λ, a vector of constraint force intensities.This leads to the following system: In the method presented in this paper, the dynamics of the air are simplified.We approximate them by considering the air pressure as uniformely distributed among the surface elements.The pressure variation inside the suction cup is imposed in our model.This is equivalent to considering an instantaneous regulation of the pressure under the cavity.This assumption is motivated by the fact that the regulation time of industrial material is usually of a few ms, which corresponds to one or two time step(s) in the simulation.
However, we do not simplify the conditions under which pressure in the cavity can be applied.In particular, at each step, our algorithm verifies that the suction cup cavity is airtight.Our pipeline proceeds as follows: 1) At the beginning of each time-step of our simulation, we compute the free-motion configuration using a preconditionned Conjugate Gradient solver.It means that we solve a first equilibrium given the internal forces and applied external forces onto the suction cup while cancelling the contact constraints.We solve the equation ( 2) with λ = 0 to obtain dv free = A −1 b.After the free motion, one can compute the violation δ free of any constraint by using its value at the beginning of the step δ 0 and the projection of the free motion into the constraint space hH(v 0 + dv free ), leading to δ free = δ 0 + hH(v 0 + dv free ).2) Then, we classify the points on the inner surface of the suction cup in case of contact.Our algorithm uses the collision response coming from the narrow-phase collision detection as input.The classification is useful to detect whether the cavity is airtight or not, in order to decide if a pressure constraint has to be solved in addition to the contact constraints.This algorithm is called the classifier in the remainder of the paper (see Section III-B).3) Once the classification is finished, we build the constraint matrix H using both the classification data and the collision response (see Section III-C).4) Finally, we compute the Schur complement of equation ( 2) W = HA −1 H T containing the coupling between the different contact points but also between the contact points and the pressure constraint.Thanks to W we are able to solve the constraints by finding their respective Lagrange multipliers.To this end, we use a Gauss-Seidel solver with successive overrelaxation (SOR) (see Section III-C).5) At the end, we apply the resulting correction forces H T λ on the points having constraints, thus applying a motion correction on the suction cup.The outcome is a new material configuration without any interpenetration while the depressurization is taken into account.

B. Classification of Constrained Zones
The suction phenomenon occurs when there is a difference between the pressure inside the cavity of the suction cup and the external pressure.It first requires the cavity of the suction cup to be airtight.Our classifier has the ability to detect this by classifying the points of the discretized inner surface of the suction cup.After computation of the free-motion configuration, the narrow-phase of the collision detection provides the signed interpenetration distances of the contact points.These distances correspond to the contact constraint violations δ free .The collision response data are then used by our classifier for distinguishing 5 classes: As explained in Figure 2, our classification works by following three consecutive stages: 1

C. Constraint Formulation and Resolution
At each time step, once the free-motion configuration of the suction cup as well as the classification have been computed, the next stage is to compute the constraint matrix H.In terms of structure, H is a sparse matrix and its embedded elements are 3D vectors (see Figure 3).H contains all the contact constraint directions as well as a field of 3D vectors which is related to the internal air pressure we choose to model as a constraint.
Each point located on the inner surface and classified as STRICTLY INSIDE and INNER BORDER is constrained by the vaccum negative pressure constraint.Given the fact that the pressure is considered as uniform inside the cavity created by the suction cup, the unique Lagrange multiplier λ p of the pressure constraint corresponds to the scalar pressure.Then H T p λ p represents the forces applied to the points which belong to the cavity.Thus, H p is a line matrix carrying inner cavity surface normals directions which are scaled by the local area of their surrounding surface elements.Assuming the surface elements are triangles, the scale factor s for an individual pressure constraint direction is evaluated as s = The constraint problem contains complementarity constraints (contact) and an equality constraint (air pressure), this MCP (Mixed Complementarity Problem) is solved with a slightly modified version of the Gauss-Seidel (GS) Algorithm described in [14].This algorithm computes the value of a vector λ so that the violations of the constraints δ = δ free + Wλ vanish.With the modified GS solver, the MCP is computed in three steps: 1) An initial guess for the Lagrange Multipliers of the contacts points is computed through the GS adapted for contact constraints, while ignoring the pressure constraint; 2) Knowing the regulated pressure of vaccum p, and the initial value of λ, we progressively increment an internal pressure λ p until it reaches p, and we reclassify the points at the same time while updating the constraint matrix H.The outcome is a fresh classification which takes the pressure into account; 3) From the updated classification and given λ p = p, the Lagrange Multipliers of the contact points which are still in contact are updated by starting again the GS solver.

D. Implementation and Performances
Our model was implemented in C++ with SOFA, an opensource simulation framework for deformable objects [15].We ran the simulations on a laptop with an Intel Core i7-7820HQ CPU at 2.90GHz, 32GB of DDR4 system memory, and a NVIDIA Quadro P3000 Mobile graphics card.We use an asynchronous preconditioner [16] with a GPU implementation in order to accelerate the construction of the W matrix.
The Conjugate Gradient solver computes the free-motion configuration using a single iteration thanks to its preconditionner, whereas we fixed a tolerance of 10 − 9 with a limit of 100 iterations for the MCP solver (Gauss Seidel).Knowing that the tetrahedral mesh of our suction cup is composed of 704 finite elements, we reach an average of 71.3 FPS when there is no contact, and we obtain an average of 11.7 FPS during suction contact.This performance drop is due to the numerous constraints which have to be solved.At a time-step inducing 274 contact constraints plus the pressure constraint inside the MCP, the measured computation costs are: 0.2 ms for the classification, 1.1 ms for the free motion computation, 14.4 ms for the construction of the W matrix, and 64.0 ms for solving the MCP.

IV. VALIDATION OF THE MODEL
Besides the performance of the simulations, we performed an evaluation based on a comparison with real data.We first validated the geometry of the suction cups that we obtained after deformation with our model.Then, we also validated the forces computed from our simulations.

A. Real Data and Experimental Setup
We designed two different types of silicone suction cups.We cast 3D printed suction cups, using respectively Dragon Skin A20 silicon ( Young modulus = 338 kPa) and Dragon Skin A30 silicon (Young modulus = 593 kPa).In the following, we will refer to the names A20 and A30 for the two suction cups.We generated the meshes using the CAD models built for 3D printing (see Figure 4).In our experimental setup, the suction surface is positionned on a flat surface, on top of a small hole which is connected to the vaccum via a tube.The pressure inside the suction cup is controlled using a pressure regulator and an associated numeric pressure sensor along the pipe tube.In this way, we were able to control the shape of the cup by gently varying the internal pressure through the regulator.

B. Geometry validation
The objective of the first experiment was to compare the shape of the real suction cup with the shape of the virtual one, in order to validate the geometry under a static approach.We applied an internal depressurization of −11 kPa and −6kPa respectively to the suction cups A30 and A20 and we recorded the initial shapes of the real suction cups as well as the final ones.
We extracted the shapes of the real suction cups using a photogrammetry procedure.For that purpose, we used the Meshroom software (https://alicevision. github.io/) to reconstruct the 3D shapes from pictures we took during the experiment (see Figure 5).Once the 3D real data were reconstructed, we could compare them with the resulting shapes from our simulation.
In order to avoid the noise resulting from the photogrammetry reconstruction, we applied smoothing and filtering operations afterwards.As illustrated in Figure 6, we compared the real suction cup with the virtual one using the Hausdorff distance.When comparing the final states, we obtained the following measurements: for the suction cup A20, we measured a Hausdorff distance of 1.59 mm.For the suction cup A30, the Hausdorff distance was 3.29 mm.

C. Force validation
The objective of the second experiment was to measure the intensity of the pull force required to detach the suction cup from the surface, assuming it was completely stuck.
As illustrated in Figure 7, we attached the top of the suction cup to an actuator with a hook.In addition, a pressure sensor was positionned between the tip of the actuator and the hook.In this setup, we consequently were able to measure the pulling force.The actuator had only one degree of freedom which means it could exclusively move vertically.As for the first experiment, we connected the cavity of the suction cup to the vaccum and were thus able to regulate the internal pressure.
We ran the experiments with both suction cups A20 and A30.For each one, we proceeded as follows: • Thanks to the actuator, the suction cup was positionned on top of the hole without pushing it too much against the flat surface.
• The vaccum was switched on and we set the cavity depressurization to a fixed value using the pressure regulator, achieving a deformation of the suction cup.• We pulled the suction cup by moving the actuator up at a speed of 5 mm/s.We reproduced a similar scenario with our simulations.We were then able to compare the force variation as a function of time and to detect the necessary pull force to detach the suction cup from the flat surface.The forces comparison between the virtual one and the ground truth to reach detachment is provided in Figure 8 with suction cup A20 and different internal pressures.

V. CONCLUSION
In this paper, we proposed a novel interactive physicallybased approach for simulating active suction cup phenomena using a constraint-based formulation.The suction cup is simulated using the FEM while the pressure is modeled with a constraint under the suction cup.We also presented an algorithmic implementation of our approach which handles the classification of the contact when a pressure variation occurs under the suction cup.Our simulation can be performed at interactive time while handling the contact points due to the suction phenonemon.We validate the performance of our approach against real data, evaluating its ability to simulate the shape as well as the forces applied to a real suction cup.As perspectives, we would like to test for different shapes for the suction cups, thus increasing the panel of suction cups that we are able to simulate with our approach.We would also like to conduct further experiments for evaluating the simulation under various boundaries conditions such as non flat surfaces or other types of contacts with the environments.Another interesting experiment would be to add friction constraints into our model in order to measure the impact of friction with regards to the suction cup detachment.Taken together, our results as well as the perspectives pave the way for improving the design as well as the control of real suction cups, especially when used in a robotics context.

Fig. 1 .
Fig. 1.Illustration of our constraint-based physically-based approach for simulating active suction cup phenomenon: (Top) the suction cup is actively stuck to the surface, (Bottom) is is then release until being completely in the air.
The first stage separates the points which are IN ACTIVE CONTACT and those which are not.It also makes an initial guess of the OUTER BORDER.2) The second stage moves the OUTER BORDER towards the region of points which are IN ACTIVE CONTACT.At the same time, it fills the region of points which are STRICTLY OUTSIDE.3) The third stage looks for an unclassified region of points and fills it as a STRICTLY INSIDE region.It also detects the associated INNER BORDER.It assumes the object has only one cavity.In addition to detecting whether the cavity of the suction cup is airtight, the classifier provides the points where the depressurization applies.We use these output data to build the constraint matrix H (Section III-C).More details about the classifier can be found in Algorithm 1.

Fig. 2 .
Fig. 2. The three consecutives stages of the classification algorithm, from top to bottom: (a) contact points are detected (red) while having an initial guess of the outer border (magenta); (b) The outer border moves towards active contact points while tagging in yellow some points that are not active on its way.The set of magenta points form the outer border; (c) The same principle is applied to detect the inner border (blue) and the set of points inside the cavity (green).

e=n e=0 ae 3
with a e the area of a triangle element e among the n surrounding elements.Each point located on the contact surface classified as IN ACTIVE CONTACT, OUTER BORDER and IN-NER BORDER creates one constraint along the normal direction of contact.When gathering the contribution of all contacts on all the nodes of the FEM, we obtain H T c λ c .Consequently, H c contains as much lines as the number of contact points.The matrices H = H c H p and H T = H T c H T p are the gathering of the two matrices H c and H p .

11 b ← OUTER BORDER 12 /
input : E contains the element pairs coming from the response of the narrow-phase collision detection output: C contains the classification of the points at the end of the algorithm execution 1 func ClassifyFromCollisionResponse(E): 2 // --STAGE 1 --3 foreach c ∈ C do 4 c ← NOT CLASSIFIED 5 foreach e ∈ E do 6 Get c ∈ C the point corresponding to e 7 c ← IN ACTIVE CONTACT 8 foreach b ∈ B do 9 if b = IN ACTIVE CONTACT then 10 stack.push(b)/ --STAGE 2 --13 while stack is not empty do 14 Let c from stack.pop() 15 c ← STRICTLY OUTSIDE 16 Get K ⊆ C the neighboor points of c is not empty do 33 Let c from stack.pop() 34 Get K ⊆ C the neighboor points of c 35 foreach k ∈ K do 36 if k = IN ACTIVE CONTACT then 37 k ← INNER BORDER 38 return C Algorithm 1: Classification of the cavity points using narrow-phase collision response data as an input

Fig. 3 .
Fig. 3.When a collision occurs, two types of 3D vectors are embedded in our constraint matrix H.The red ones are the contact constraint directions whereas the green ones constitute the uniform negative pressure distribution in the cavity.The norms of the latters are proportionnal to their respective surrounding surface elements.

Fig. 6 .
Fig.6.Shape comparaison between the real suction cup (grey) and the virtual one (red) at their final states.

Fig. 7 .
Fig. 7. Experimental setup for the force measurements.The suction cup is attached to a force sensor.When it is positionned on a flat surface, its cavity is linked to a vaccum pump with a regulator inbetween.

Fig. 8 .
Fig. 8. Force measurement from real data and simulations (called virtual) for different depressurizations.The force increases because the actuator pulls the suction to the top.The force drops to zero when the suction cup distachs from a flat surface.Simulation forces are represented with dashed lines.