An effective nonlinear dynamic formulation to analyze grasping capability of soft pneumatic robotic gripper

ABSTRACT

Soft pneumatic robotic grippers have found extensive applications across various engineering domains, which prompts active research due to their splendid compliance, high flexibility, and safe human-robot interaction over conventional stiff counterparts. Previously simplified rod-based models principally focused on clarifying overall large deformation and bending postures of soft grippers from static or quasi-static perspectives, whereas it is challenging to elaborate grasping characteristics of soft grippers without considering contact interaction and nonlinear large deformation behaviors. To address this, based on absolute nodal coordinate formulation (ANCF), comprehensively allowing for structural complexity, geometric, material and boundary nonlinearities, and incorporating Coulomb’ friction law with a multiple-point contact method, we put forward an effective nonlinear dynamic modeling approach for exploring grasping capability of soft gripper. Moreover, we solved the established dynamic equations using Generalized-α scheme, and conducted thorough numerical simulation analysis on a three-jaw soft pneumatic gripper (SPG) in terms of grasping configurations, displacements and contact forces. The proposed dynamic approach can accurately both describe complicated deformed configurations along with stress distribution and provide a feasible solution to simulate grasping targets, whose effectiveness and precision were analyzed theoretically and verified experimentally, which may shed new light on devising and optimizing other multifunctional SPGs.

GRAPHICAL ABSTRACT

The configuration variations and von Mises stress distribution of the SPG at different moments.

KEYWORDS:

• Nonlinear dynamic formulation
• soft pneumatic robotic grippers
• contact interaction
• grasping capability
• configuration variations

1. Introduction

The advancement of soft materials such as silicone, hydrogel, electroactive polymer, shape memory alloy and ionic polymer metal composite [1–3], has significantly accelerated the development of diverse soft grippers which have huge application value in all walks of life. Over the past two decades, soft pneumatic grippers (SPGs) are increasingly gaining attention owing to their merits of excellent flexibility, large deformation, strong robustness, high safety and good environmental compatibility [4]. The design, materials and control of SPGs generally follow a certain basic standard, and completely compliant soft grippers composed of embedded pneumatic networks in soft materials are usually actuated by compressed air. Application of air pressures to chambers of soft grippers induces shape changes accordingly and generates diverse deformation styles including stretching, contracting, bending, twisting, curling, twining, and combinations thereof [5, 6]. Consequently, SPGs can achieve not only simple 2D motion but also complex 3D motion, enabling dexterous actions for grasping flexible and fragile objects with variable sizes or irregular shapes in a non-damaging manner.

Currently, a wide variety of investigations have primarily concentrated on the design, fabrication, and grasping experiments of different SPGs tailored for some specific tasks. Representative examples include a starfish-like SPG [7], an octopus-shaped tentacle [8], a five-finger anthropomorphic SPG [9], a four-finger SPG [10–12], a three-finger SPG [13], and other SPGs [14–18]. However, how to model and simulate grasping motion of these SPGs remain a topic worthy of exploration.

To date, numerous research works have addressed static or kinematic modeling and simulation analysis of SPGs. Zhou et al. [19] analyzed three different point contact grasping modes of soft hands, and constructed a static model to study the intrinsic curvature profile and grasping stability. Payrebrune and O’Reilly [20] developed a rod-based model with five-parameter constitutive relation for soft robotic arms to predict their bending deformation postures, and explored intrinsic curvature-pressure and moment-curvature relations through a comparison of finite element model. Uppalapati and Krishnan [21] proposed a pure helical model and a spatial Cosserat rod model for a spiral SPG to capture its deformed behaviors and motion characteristics. Hao et al. [22] built a mathematical model that captures the relationship between geometrical properties of chambers and applied pressures to depict bending deformation of a soft gripper, and analyzed the bidirectional bending extent through quasi-static nonlinear simulation. Gong et al. [23] established a kinematic model for a soft robotic arm with 3D movement based on homogeneous transformation matrices, and predicted its location and motion through mathematical modeling. Wu et al. [24] developed a static model, performed kinematic analysis on a SPG with a jointed endoskeleton structure and revealed a quasi-linear relation between grasping force and input pressure. Based on minimum potential energy method and continuum rod theory, Gu et al. [25] exploited an analytical model for a bionic helical SPG to inspect the influences of geometric parameters, material properties and external force on 3D deformations. Lotfiani et al. [26] presented an analytical model to predict bending, contact force, twist angle of a SPG, and validated the accuracy of theoretical model through finite element simulations and experiments. Besides, there are also lots of relevant investigations on modeling and motion simulation of SPGs [27–29], which will not be listed out here.

To sum up, most previous investigations on SPGs almost focused on design, fabrication, quasi-static or static modeling, simulation analysis, and grasping experiments. However, some studies approximately treated soft grippers as conventional multi-segment rigid rods, and employed Lagrangian method for dynamic modeling and analysis, but it is difficult to establish accurate models and reveal grasping performance and nonlinear dynamic response [24, 30, 31]. Moreover, there is a scarcity of research on dynamic modeling of SPGs considering contact interaction and nonlinear large deformation. There are also limited in-depth discussions on which methods to be adopted to build an effective dynamic model for SPGs, which is a challenging problem that needs to be addressed urgently. Building on our prior research [32], where we proposed a multiple-point contact method based on absolute nodal coordinate formulation (ANCF) for analyzing a soft pneumatic actuator with nonlinear large deformation, and validated its effectiveness and accuracy. This paper extends our investigation to propose a nonlinear dynamic modeling approach for a specific three-jaw SPG and study its grasping capabilities.

In this paper, within the framework of ANCF, uniting a multiple-point contact method with Coulomb’s friction law, an accurate dynamic model is established to analyze the grasping characteristics of the three-jaw SPG. This paper is organized as follows: Section 2 sets forth overall structural design of the SPG. Section 3 deduces different nonlinear force vectors and their derivatives, and builds the corresponding dynamic model. Section 4 implements dynamic simulation analysis and experimental investigation on grasping process of the SPG, clarifies its nonlinear deformation, configuration variations, stress distributions, and then compares theoretical results with experimental ones to verify the validity and precision of the proposed dynamic formulation. Section 5 draws several conclusions and looks ahead to potential applications. This work aims to contribute to a deeper understanding of SPGs’ grasping behavior and facilitate their application in practical engineering.

2. The structural design of a three-jaw SPG

In daily life, numerous objects, such as balloons, lifebuoy, rubber boats, and trampolines, are composed of inflatable and deformable soft materials [33]. Leveraging the principle of inflatable deformation, combining a stable tripod structure with pneumatic networks actuators, we devise a three-jaw soft pneumatic gripper (SPG) which comprises of one base and three soft pneumatic actuators (embedded pneumatic networks in silicone elastomers) positioned at 120° to each other, and three actuators are inserted into three symmetrical rectangular holes in the base according to positional pairing relationship. The bending deformation extent of the SPG in 3D space can be attained by adjusting the magnitude of the air pressure. Figures 1 (a)-(e) illustrates the detailed construction of the SPG. Referring to literature [34], and the detailed fabrication process is as follows:

Figure 1. The structural design of a three-jaw SPG. (a) a soft pneumatic bending actuator after inflation, (b) sectional view of the actuator, (c) a base, (d) the SPG and a ball in a still state, (e) the construction and components of the SPG, (f) the SPG and the ball in grasping state.

Step 1: Use Dragon Skin 30 silicone material, mix the two components A and B in a 1:1 volume ratio, and then quickly and evenly stir them in a blender;

Step 2: Apply vaseline evenly on the surface of the upper mold to facilitate subsequent demolding;

Step 3: Evenly inject the prepared silicone material into the upper mold, after the silicone is fully precipitated and filled into the upper mold, remove excess silicone, and use a vacuum pump to vacuum the upper mold filled with silicone for 5 minutes to remove air bubbles;

Step 4: Dry the upper mold filled with silicone at room temperature for 12 hours, then slowly remove the upper mold to obtain the upper air chamber;

Step 5: Drill a hole on the left end face of the upper air chamber, insert the PVC hose into the hole, and seal the interface with silicone adhesive;

Step 6: Evenly inject the prepared silicone material into the lower mold, after vacuuming for 5 minutes, place the upper air chamber on the bottom mold filled with silicone, dry at room temperature for 6 hours, and remove the mold to obtain a soft pneumatic actuator;

Step 7: Check the air tightness, put the soft pneumatic actuator into the water, use the air pump to inflate it slowly, observe whether there are air bubbles escaping;

Step 8: Insert three identical soft pneumatic actuators into a 3D printed base according to corresponding positional relationship to obtain a three-jaw soft pneumatic gripper.

Table 1. Description of specific parameters of the three-jaw SPG.

Detailed explanation	Parameters	Magnitude
Length of soft actuator	L	44 mm
Width of soft actuator	W	20 mm
Height of soft actuator	H	17 mm
Height of internal air chambers	H1	13 mm
Thickness of air chamber walls	t	1 mm
Length of solid section	t1	10 mm
Interval between adjacent contact walls	t2	2 mm
Length of each air chamber	t3	6 mm
Length of front air chamber	t4	8 mm
Side length of base	l	42 mm
Width of base	w	30 mm
Thickness of base	h	4 mm
External length of rectangular slot	L1	23 mm
Internal length of rectangular slot	L2	20 mm
External width of rectangular slot	W1	20 mm
Internal width of rectangular slot	W2	17 mm
External height of rectangular slot	h1	10 mm
Internal height of rectangular slot	h2	8 mm

When air pressure is applied to air chambers of the soft pneumatic actuator through an air passage, expansion deformation occurs. Upon reaching a certain extent, adjacent chambers come into contact to deform the soft actuator. The constraining layer restricts the axial elongation of the actuator’s plane side, and the axial elongation of air chambers becomes a dominant factor which causes the actuator to bend inward. Synchronously applying air pressure to three soft actuators enables the SPG to perform bending actions, resembling human’s fingers for grasping objects of various shapes (e.g. small ball, glass, sea cucumber). The SPG is naturally in a still state, and closed state under positive pressure. This paper mainly focuses on nonlinear grasping characteristics of the SPG and contact interaction between the SPG and target objects. Table 1 lists specific parameters of the three-jaw SPG in Figure 1.

3. Deduction of different force vectors and dynamic modeling

3.1. The element elastic force vector and its derivative

While the previously proposed simplified rod-based model is generally useful for describing global bending deformation and postures of soft grippers, it falls short in revealing nonlinear contact characteristics and grasping behaviors between grippers and target objects. To address this limitation, we present a dynamic modeling approach by combining the hyperelastic constitutive Yeoh model with ANCF solid elements [35]. This approach can be extended to explore grasping mechanisms of complex soft grippers. Next, we will describe nonlinear modeling process and derive several different force vectors detailedly.

Based on a general theory of flexible bodies [36, 37], uniting eight-node ANCF hexahedron solid elements with the selective reduced integration technique to relieve volumetric locking problem of soft bodies [38]. Under this circumstances, ANCF solid elements can accurately capture soft grippers' large bending deformation.

The absolute position coordinate vector of an arbitrary point in an ANCF solid element in Cartesian coordinate system can be expressed as

(1) $\mathbf{r}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta },\mathit{t}\right)=\mathbf{S}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta }\right){\mathbf{q}}_e(\mathit{t})$(1)

where $\mathit{S}=\left[\begin{array}{cc}{S}_1{\mathbf{I}}_{\mathbf{3}}& S_2{\mathbf{I}}_{\mathbf{3}}\end{array}\begin{array}{cc}\dots & S_8{\mathbf{I}}_{\mathbf{3}}\end{array}\right]$ is an element shape function matrix, I₃ is a 3 × 3 identity matrix, $(\mathit{\xi },\mathit{\eta },\mathit{\zeta })$ is a parameterized coordinate vector defined in a natural coordinate system, and $-1\le \mathit{\xi }\le 1$, $-1\le \mathit{\eta }\le 1$, $-1\le \mathit{\zeta }\le 1$. ${\mathbf{q}}_e(\mathit{t})$ is time-dependent nodal position coordinate vector of the ANCF solid element. Based on a generalized form of the shape functions for the solid element, the specific shape function is formulated as

(2) $S_i(\mathit{\xi },\mathit{\eta },\mathit{\zeta })=\frac{1}{8}\left(1+{\xi }_{\mathrm{i}}\mathit{\xi }\right)\left(1+{\eta }_{\mathrm{i}}\mathit{\eta }\right)\left(1+{\zeta }_{\mathrm{i}}\mathit{\zeta }\right),\mathit{i}=1,........8$(2)

The transformation relationship between the derivatives of the shape functions with respect to Cartesian coordinate and natural coordinate is denoted as

(3) ${{\left[\begin{array}{ccc}\mathrm{\partial }{S}_i/\mathrm{\partial x}& \mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{y}& \mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{z}\end{array}\right]}^{\mathit{\hspace{0.17em}}}}^{\mathrm{T}}={\mathbf{J}}^{-1}{{\left[\begin{array}{ccc}\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\xi }& \mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\eta }& \mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\zeta }\end{array}\right]}^{\mathit{\hspace{0.17em}}}}^{\mathrm{T}}$(3)

where $\mathbf{J}$ is Jacobian matrix and written as

(4) $\mathbf{J}=\left[\begin{array}{ccc}\sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\xi }\cdot x_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\xi }\cdot y_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\xi }\cdot z_i\\ \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\eta }\cdot x_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\eta }\cdot y_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\eta }\cdot z_i\\ \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\zeta }\cdot x_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\zeta }\cdot y_i& \sum _{\mathit{i}=1}^8\mathrm{\partial }{S}_i/\mathrm{\partial }\mathit{\zeta }\cdot z_i\end{array}\right]$(4)

By introducing the volumetric energy penalty function [39], the total strain energy of Yeoh model for hyperelastic incompressible soft elastomers is deduced as

(5) $U_Y={\int }_V\left[\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}{\left(I_1-3\right)}^{\mathit{i}}\right]\mathit{dV}+\frac{1}{2}\mathit{\kappa }{\left(J_C-1\right)}^2\mathit{V}$(5)

Here, $J_C$ is the determinant of the deformation gradient tensor ${\mathbf{J}}_C=\left[\begin{array}{ccc}\mathrm{\partial }{\mathbf{r}}_C/\mathrm{\partial }\mathit{x}& \mathrm{\partial }{\mathbf{r}}_C/\mathrm{\partial }\mathit{y}& \mathrm{\partial }{\mathbf{r}}_C/\mathrm{\partial }\mathit{z}\end{array}\right]$ at the center point C of the ANCF solid element, which is signified as $J_C=\det ({\mathbf{J}}_C)$. $I_1$ is the trace of the right Cauchy-Green deformation tensor $\mathbf{C}={\mathbf{J}}_C^{\mathrm{T}}{\mathbf{J}}_C$, and $\mathit{\kappa }$ is an incompressibility constant which shall be selected adequately large to conform with incompressibility condition $J_C=1$. Three material constants C_i0 (i = 1, 2, 3) can be determined by uniaxial tensile test [39].

Differentiating Equation (5)(5) $U_Y={\int }_V\left[\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}{\left(I_1-3\right)}^{\mathit{i}}\right]\mathit{dV}+\frac{1}{2}\mathit{\kappa }{\left(J_C-1\right)}^2\mathit{V}$(5) relative to the element position coordinate vector ${\mathbf{q}}_e$, an element elastic force vector is given by ${\mathbf{Q}}_{\mathit{ek}}\mathit{\hspace{0.17em}}=\mathrm{\partial }{U}_Y/\mathrm{\partial }{\mathbf{q}}_e$ and derived as

(6) ${\mathbf{Q}}_{\mathit{ek}}={\int }_V\left[\sum _{\mathit{i}=1}^3\mathit{i}{C}_{\mathit{i}0}{\left(I_1-3\right)}^{\mathit{i}-1}\right]{\left(\frac{\mathrm{\partial }{I}_1}{\mathrm{\partial }{\mathbf{q}}_e}\right)}^{\mathrm{T}}\mathit{dV}+\mathit{\kappa }\mathit{V}\left(J_C-1\right)\mathit{\hspace{0.17em}}{\left(\frac{\mathrm{\partial }{\mathit{J}}_{\mathrm{C}}}{\mathrm{\partial }{\mathbf{q}}_e}\right)}^{\mathrm{T}}$(6)

Differentiating the element elastic force vector ${\mathbf{Q}}_{\mathit{ek}}\mathit{\hspace{0.17em}}$ relative to the element absolute coordinate vector ${\mathbf{q}}_e$, and thus the corresponding derivative is expressed as

(7) $\begin{array}{l}{K}_{\mathit{ek}}=\frac{\partial {\mathit{Q}}_{\mathit{ek}}}{\partial q_e}=\\ {\displaystyle {\int }_V\left\{\begin{array}{l}{C}_{10}\frac{\partial }{\partial q_e}{\left(\frac{\partial I_1}{\partial q_e}\right)}^{\text{T}}+2C_{20}\left(I_1-3\right)\frac{\partial }{\partial q_e}{\left(\frac{\partial I_1}{\partial q_e}\right)}^{\text{T}}+2C_{20}{\left(\frac{\partial I_1}{\partial q_e}\right)}^{\text{T}}\frac{\partial I_1}{\partial q_e}\\ +3C_{30}{\left(I_1-3\right)}^2\frac{\partial }{\partial q_e}{\left(\frac{\partial I_1}{\partial q_e}\right)}^{\text{T}}+6C_{30}\left(I_1-3\right){\left(\frac{\partial I_1}{\partial q_e}\right)}^{\text{T}}\frac{\partial I_1}{\partial q_e}\end{array}\right\}}\mathit{dV}\\ +\mathit{\kappa V}\left(J_C-1\right)\frac{\partial }{\partial q_e}{\left(\frac{\partial J_{\text{C}}}{\partial q_e}\right)}^{\text{T}}+\mathit{\kappa V}{\left(\frac{\partial J_{\text{C}}}{\partial q_e}\right)}^{\text{T}}\left(\frac{\partial J_{\text{C}}}{\partial q_e}\right)\end{array}$(7)

3.2. The element pressure force vector and its derivative

In this paper, we assume that the air chamber is a cuboid, Figure 2 displays the surface element of air chamber walls.

Figure 2. The surface element of air chamber walls.

The virtual work done by an element pressure force is indicated as

(8) $\mathit{\delta }{W}_{\mathit{ep}}={\int }_A\mathit{\delta }{\bar{\mathbf{r}}}^{\mathrm{T}}\mathit{p}\mathbf{n}\mathit{dA}=\mathit{\delta }\bar{\mathbf{q}}{ }_e{{ }^p}^{\mathrm{T}}{\mathbf{Q}}_{\mathit{ep}}$(8)

Here, ${\mathbf{Q}}_{\mathit{ep}}$ is an element pressure force vector and represented as

(9) ${\mathbf{Q}}_{\mathit{ep}}={\int }_A{\bar{\mathbf{S}}}^{\mathrm{T}}\mathit{p}\mathbf{n}\mathit{dA}$(9)

Here, p is the value of air pressure, $\bar{\mathbf{r}}$ is a position coordinate vector of an arbitrary point on the air chamber wall in the global coordinate system, which is given by $\bar{\mathbf{r}}=\bar{\mathbf{S}}{\bar{\mathbf{q}}}_{\mathit{e}}^{\mathit{p}}$, and $\bar{\mathbf{S}}=\left[{\bar{S}}_1{\mathbf{I}}_2{\bar{S}}_2{\mathbf{I}}_2{\bar{S}}_3{\mathbf{I}}_2{\bar{S}}_4{\mathbf{I}}_2\right]$ is the shape function matrix of a rectangular surface element subjected to air pressure. In this context, I₂ is a 2 × 2 identity matrix. Moreover, ${\bar{\mathbf{q}}}_{\mathit{e}}^{\mathit{p}}$, $\mathbf{n}$, A are severally an element position coordinate vector, a unit normal vector, a surface area of the rectangular surface element, and $\mathbf{n}={\mathbf{n}}_n/||{\mathbf{n}}_n||$ is the normal vector of the deformed curved surface.

Differentiating the element pressure force vector ${\mathbf{Q}}_{\mathit{ep}}$ relative to the element position coordinate vector ${\bar{\mathbf{q}}}_{\mathit{e}}^{\mathit{p}}$, and the derivative of ${\mathbf{Q}}_{\mathit{ep}}$ is denoted as

(10) ${\mathbf{K}}_{\mathit{ep}}=\frac{\mathrm{\partial }{\mathbf{Q}}_{\mathit{ep}}}{\mathrm{\partial }\bar{\mathbf{q}}\mathit{ep}}={\int }_A{\bar{\mathbf{S}}}^{\mathrm{T}}\mathit{p}\frac{\mathrm{\partial }\mathbf{n}}{\mathrm{\partial }\bar{q}\mathit{ep}}\mathit{dA}$(10)

3.3. The element gravitational and concentrated vectors and their derivatives

Similarly, the virtual work done by the element gravitational force is written as

(11) $\mathit{\delta }{W}_{\mathit{eg}}={\int }_V\mathit{\delta }{\mathbf{r}}^{\mathrm{T}}\mathit{\rho }\mathbf{g}\mathit{dV}=\mathit{\delta }{{\mathbf{q}}_e}^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eg}}$(11)

Substituting Equation (1)(1) $\mathbf{r}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta },\mathit{t}\right)=\mathbf{S}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta }\right){\mathbf{q}}_e(\mathit{t})$(1) into Equation (11)(11) $\mathit{\delta }{W}_{\mathit{eg}}={\int }_V\mathit{\delta }{\mathbf{r}}^{\mathrm{T}}\mathit{\rho }\mathbf{g}\mathit{dV}=\mathit{\delta }{{\mathbf{q}}_e}^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eg}}$(11) , the element gravitational force vector has the following expression

(12) ${\mathbf{Q}}_{\mathit{eg}}={\int }_V{\mathbf{S}}^{\mathrm{T}}\mathit{\rho }\mathbf{g}\mathit{dV}$(12)

Then, the virtual work done by the element concentrated force is given by

(13) $\mathit{\delta }{W}_{\mathit{eF}}=\underset{\mathit{V}}{\int }\mathit{\delta }{\mathbf{r}}^{\mathrm{T}}\mathit{\rho adV}=\mathit{\delta }{\mathbf{q}}_e^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eF}}$(13)

Substituting Equation (1)(1) $\mathbf{r}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta },\mathit{t}\right)=\mathbf{S}\left(\mathit{\xi },\mathit{\eta },\mathit{\zeta }\right){\mathbf{q}}_e(\mathit{t})$(1) into Equation (13)(13) $\mathit{\delta }{W}_{\mathit{eF}}=\underset{\mathit{V}}{\int }\mathit{\delta }{\mathbf{r}}^{\mathrm{T}}\mathit{\rho adV}=\mathit{\delta }{\mathbf{q}}_e^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eF}}$(13) , the element concentrated force vector takes the following form

(14) ${\mathbf{Q}}_{\mathit{eF}}=\underset{\mathit{V}}{\int }{\mathbf{S}}^{\mathrm{T}}\mathit{\rho adV}$(14)

Here, g and $\mathbf{a}$ are the gravitational acceleration vector and the motion acceleration vector of the SPG, respectively. Since the element gravitational force vector ${\mathbf{Q}}_{\mathit{eg}}$ and the element concentrated force ${\mathbf{Q}}_{\mathit{eF}}$ are constant, the derivatives ${\mathbf{K}}_{\mathit{eg}}$, ${\mathbf{K}}_{\mathit{eF}}$ of ${\mathbf{Q}}_{\mathit{eg}}$, ${\mathbf{Q}}_{\mathit{eF}}$ relative to the element absolute coordinate vector ${\mathbf{q}}_e$ are both equal to 0. Thus, the derivatives ${\mathbf{K}}_g$ and ${\mathbf{K}}_F$ after assembling also equal to 0.

3.4. The element contact force vector and its derivative

The previously developed multiple-point contact method, in which multiple Gaussian points of a slave element contact with projection points of a master element, is adopted to gain contact forces [32]. The contact zone is approximated as a group of integration points which project onto the corresponding master element. Figure 3 displays the closest point projection in three-dimensional coordinate system. s is a slave point on a slave surface, and m is a master point on a master surface which is the projection point from the slave point s onto the master surface [40]. ${\mathbf{r}}_s$ is the global position vector of the slave point s, ${\mathbf{r}}_m$ is the global position vector of the master point m. Defining ${\mathbf{r}}_{\mathit{ms}}={\mathbf{r}}_s-{\mathbf{r}}_m$ as the relative position vector of the slave point s relative to the master point m. The closest point projection procedure [41] induces a following minimum problem.

(15) $\parallel {\mathbf{r}}_{\mathit{ms}}\parallel ={\mathbf{r}}_{ms}^{\mathrm{T}}{\mathbf{r}}_{\mathit{ms}}\to \min$(15)

Figure 3. Closet point projection procedure.

$g_n$ is defined as the dot product of the relative position vector ${\mathbf{r}}_{\mathit{ms}}$ and the corresponding normal vector ${\mathbf{n}}_m$, that is, $g_n={\mathbf{n}}_m^{\mathrm{T}}{\mathbf{r}}_{\mathit{ms}}$. It is obvious that for $g_n>0$, the slave point s does not penetrate into the master surface, meaning there is no contact. As for $g_n\le 0$, it is determined that the slave point s penetrates into the master surface, indicating that penetration has occurred. The value of penetration is ${\delta }_n=\left|g_n\right|=-g_n$, and the normal contact force per unit area is formulated as

(16) $p_n=k_n{\delta }_n$(16)

where $k_n$ is the normal contact stiffness per unit area. Base on the penalty function method, the normal contact tractions acting on an arbitrary point of the slave surface and the corresponding projection point of master surface take the form

(17) ${\mathbf{p}}_s=k_n{\delta }_n{\mathbf{n}}_m,{\mathbf{p}}_m=-k_n{\delta }_n{\mathbf{n}}_m$(17)

The position variation of the slave point s and the master point m are given by

(18) $\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}={\bar{\mathbf{S}}}_{\mathit{s}}\mathit{\delta }\bar{\mathbf{q}}\mathit{es},\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}={\bar{\mathbf{S}}}_{\mathit{m}}\mathit{\delta }\bar{\mathbf{q}}\mathit{em}$(18)

Here, ${\bar{\mathbf{S}}}_{\mathit{s}},{\bar{\mathbf{S}}}_{\mathit{m}}$ are severally the shape function matrices of the slave and master elements of the contact pair, and $\bar{\mathbf{q}}\mathit{es},\bar{\mathbf{q}}\mathit{em}$ severally signify the element position coordinate vectors of the slave and master elements.

The virtual work done by the normal contact tractions vectors is expressed as

(19) $\mathit{\delta }{W}_{\mathit{ec}}={\int }_{\mathit{As}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}^{\mathrm{T}}{\mathbf{p}}_s\mathit{dA}+{\int }_{\mathit{Am}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}^{\mathrm{T}}{\mathbf{p}}_m\mathit{dA}$(19)

Substituting Equation (17)(17) ${\mathbf{p}}_s=k_n{\delta }_n{\mathbf{n}}_m,{\mathbf{p}}_m=-k_n{\delta }_n{\mathbf{n}}_m$(17) and Equation (18)(18) $\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}={\bar{\mathbf{S}}}_{\mathit{s}}\mathit{\delta }\bar{\mathbf{q}}\mathit{es},\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}={\bar{\mathbf{S}}}_{\mathit{m}}\mathit{\delta }\bar{\mathbf{q}}\mathit{em}$(18) into Equation (19)(19) $\mathit{\delta }{W}_{\mathit{ec}}={\int }_{\mathit{As}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}^{\mathrm{T}}{\mathbf{p}}_s\mathit{dA}+{\int }_{\mathit{Am}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}^{\mathrm{T}}{\mathbf{p}}_m\mathit{dA}$(19) , the virtual work is rephrased as

(20) $\mathit{\delta }{W}_{\mathit{ec}}=\mathit{\delta }\bar{\mathbf{q}}{ }_e^s{ }^{\mathrm{T}}{\mathbf{Q}}_{ec}^s+\mathit{\delta }\bar{\mathbf{q}}{ }_e^m{ }^{\mathrm{T}}{\mathbf{Q}}_{ec}^m$(20)

where the element contact force vectors of the slave and master elements are represented as

(21) ${\mathbf{Q}}_{ec}^s={\int }_{\mathit{As}}{\bar{\mathbf{S}}}_{\mathit{s}}^{\mathrm{T}}{k}_n{\delta }_n{\mathbf{n}}_m\mathit{dA},{\mathbf{Q}}_{ec}^m=-{\int }_{\mathit{Am}}{\bar{\mathbf{S}}}_{\mathit{m}}^{\mathrm{T}}{k}_n{\delta }_n{\mathbf{n}}_m\mathit{dA}$(21)

Differentiating the element contact force vectors ${\mathbf{Q}}_{ec}^s,{\mathbf{Q}}_{ec}^m$ relative to the corresponding element position coordinate vectors $\overline{q}\mathit{es},\overline{q}\mathit{em}$, and thus the derivatives of the element contact force vectors are severally denoted as

(22) $K_{ec}^s=\frac{\partial Q_{ec}^s}{\partial \overline{q}\mathit{es}},K_{ec}^m=\frac{\partial Q_{ec}^m}{\partial \overline{q}\mathit{em}}$(22)

3.5. The element frictional force vector and its derivative

Figure 4 (a) displays the configuration of a contact pair with frictional contact. s is a slave point in a slave surface, and m is the corresponding projection point from the slave point s onto the master surface. ${\mathbf{n}}_m$ is the unit normal vector of the contact surface at the projection point m, ${\mathbf{r}}_s$ is the global position vector of the slave point s, and ${\mathbf{r}}_m$ is the global position vector of the projection point m.

Figure 4. The contact pair. (a) Configuration of contact pair with frictional contact, (b) Tangential slipping of the projection point m on the master surface.

Figure 4 (b) displays the tangential slipping of the projection point m on the master surface. The tangential frictional force can be gained by computing the tangential slipping displacement Δg_τ of the projection point m on the master surface [41]. Assume that m_t and m_T are severally the projection points on the master surface at the current moment t and the latter moment t+Δt. The tangential slipping displacement Δg_τ can be indicated as

(23) $\mathrm{\Delta }{\mathbf{g}}_{\tau }={\mathbf{r}}_{m_T}-{\mathbf{r}}_{m_t}$(23)

Define the frictional force vector per unit area acting on the slave point s as ${\mathbf{f}}_s=\mathbf{f}$. On the basis of the above derivation process, it is known from Coulomb’ friction law that if the norm $||\mathbf{f}||$ is less than or equal to the maximum static frictional force $f_{\max }={\mu }_s{k}_n{\delta }_n$, where ${\mu }_s$ is the static frictional coefficient and ${\delta }_n$ is the value of penetration, there is no relative slipping between contact interfaces in the tangential direction, and the contact pair is in the stick phase. Inversely, if the norm $||\mathbf{f}||$ is greater than the maximum static frictional force $f_{\max }$, there is relative slipping between contact interfaces, and the contact pair is in the slip phase, and the magnitude of the frictional force is proportional to the normal contact force and its direction is opposite to the direction of the relative velocity.

In the stick phase, the variation of $\mathbf{f}$ is proportional to the displacement variation $\mathrm{\Delta }{\mathbf{g}}_{\tau }$ of the projection point m from the current moment t to the latter moment t+Δt, and their directions are opposite, which is expressed as

(24) $\mathrm{\Delta }\mathbf{f}\left(\mathit{t}\right)=-k_{\tau }\mathrm{\Delta }{\mathbf{g}}_{\tau }\left(\mathit{t}\right)$(24)

where k_τ is the tangential contact stiffness per unit area.

Defining $\mathbf{f}(\mathit{t}+\mathrm{\Delta t})$ and $\mathbf{f}(\mathit{t})$ as the frictional force vectors per unit area acting on the slave point s at the latter moment t+Δt and the current moment t, respectively. The tangential component of $\mathbf{f}(\mathit{t})$ at the current moment t is ${\mathbf{f}}_{\tau }\left(\mathit{t}\right)=\mathbf{f}\left(\mathit{t}\right)-{\mathbf{f}}_n\left(\mathit{t}\right)$, ${\mathbf{f}}_n\left(\mathit{t}\right)$ is the normal component of $\mathbf{f}\left(\mathit{t}\right)$ and given by ${\mathbf{f}}_n\left(\mathit{t}\right)={\mathbf{n}}_m{\mathbf{n}}_m^{\mathrm{T}}\mathbf{f}\left(\mathit{t}\right)$, and the relationship between $\mathbf{f}\left(\mathit{t}+\mathrm{\Delta t}\right)$ and $\mathbf{f}\left(\mathit{t}\right)$ is written as

(25) $\mathbf{f}\left(\mathit{t}+\mathrm{\Delta t}\right)={\mathbf{f}}_{\tau }\left(\mathit{t}\right)+\mathrm{\Delta }\mathbf{f}\left(\mathit{t}\right)=\left(\mathbf{I}-{\mathbf{n}}_m{\mathbf{n}}_m^{\mathrm{T}}\right)\mathbf{f}\left(\mathit{t}\right)-k_{\tau }\mathrm{\Delta }{\mathbf{g}}_{\tau }\left(\mathit{t}\right)$(25)

From Equation (25)(25) $\mathbf{f}\left(\mathit{t}+\mathrm{\Delta t}\right)={\mathbf{f}}_{\tau }\left(\mathit{t}\right)+\mathrm{\Delta }\mathbf{f}\left(\mathit{t}\right)=\left(\mathbf{I}-{\mathbf{n}}_m{\mathbf{n}}_m^{\mathrm{T}}\right)\mathbf{f}\left(\mathit{t}\right)-k_{\tau }\mathrm{\Delta }{\mathbf{g}}_{\tau }\left(\mathit{t}\right)$(25) , the tangential frictional force vector per unit area acting on the slave point s is represented as

(26) ${\mathbf{f}}_s=\left\{\begin{array}{c}\mathbf{f},\parallel \mathbf{f}\parallel f_{max}\\ {\mu }_d{k}_n{\delta }_n\frac{\mathbf{f}}{\parallel \mathbf{f}\parallel },\parallel \mathbf{f}\parallel >f_{max}\end{array}\right.$(26)

where μ_d is the dynamic frictional coefficient.

Based on the relationship between action force and reaction force, and the frictional force vector per unit area acting on the projection point m is denoted as

(27) ${\mathbf{f}}_m=-{\mathbf{f}}_s$(27)

The virtual work done by the element frictional force vectors takes the form

(28) $\mathit{\delta }{W}_{\mathit{ef}}={\int }_{\mathit{As}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}^{\mathrm{T}}{\mathbf{f}}_s\mathit{dA}+{\int }_{\mathit{Am}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}^{\mathrm{T}}{\mathbf{f}}_m\mathit{dA}$(28)

According to $\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}={\bar{\mathbf{S}}}_{\mathit{s}}\mathit{\delta }\bar{\mathbf{q}}\mathit{es},\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}={\bar{\mathbf{S}}}_{\mathit{m}}\mathit{\delta }\bar{\mathbf{q}}\mathit{em}$, substituting Equation (26)(26) ${\mathbf{f}}_s=\left\{\begin{array}{c}\mathbf{f},\parallel \mathbf{f}\parallel f_{max}\\ {\mu }_d{k}_n{\delta }_n\frac{\mathbf{f}}{\parallel \mathbf{f}\parallel },\parallel \mathbf{f}\parallel >f_{max}\end{array}\right.$(26) and Equation (27)(27) ${\mathbf{f}}_m=-{\mathbf{f}}_s$(27) into Equation (28)(28) $\mathit{\delta }{W}_{\mathit{ef}}={\int }_{\mathit{As}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{s}}^{\mathrm{T}}{\mathbf{f}}_s\mathit{dA}+{\int }_{\mathit{Am}}\mathit{\delta }{\bar{\mathbf{r}}}_{\mathit{m}}^{\mathrm{T}}{\mathbf{f}}_m\mathit{dA}$(28) , and the virtual work is rephrased as

(29) $\mathit{\delta }{W}_{\mathit{ef}}=\mathit{\delta }\overline{q}\mathit{es}\text{T}{Q}_{\text{ef}}^{\text{s}}\text{+}\mathit{\delta }{\overline{q}}_{\mathit{e}}^{\mathit{m}\text{T}}{Q}_{\text{ef}}^{\text{m}}$(29)

where the element tangential frictional force vectors of the slave and master elements are formulated as

(30) ${\mathbf{Q}}_{ef}^s={\int }_{\mathit{As}}{\bar{\mathbf{S}}}_{\mathit{s}}^{\mathrm{T}}{\mathbf{f}}_s\mathit{dA},{\mathbf{Q}}_{ef}^m={\int }_{\mathit{Am}}{\bar{\mathbf{S}}}_{\mathit{m}}^{\mathrm{T}}{\mathbf{f}}_m\mathit{dA}$(30)

Differentiating the element tangential frictional force vectors ${\mathbf{Q}}_{ef}^s,{\mathbf{Q}}_{ef}^m$ relative to the element position coordinate vectors $\bar{\mathbf{q}}\mathit{es},\bar{\mathbf{q}}\mathit{em}$, and thus the derivatives of the corresponding element tangential frictional force vectors are severally indicated as

(31) $K_{ef}^s=\frac{\partial Q_{ef}^s}{\partial \overline{q}\mathit{es}},K_{ef}^m=\partial Q_{ef}^m\partial \overline{q}\mathit{em}$(31)

Based on the above derivation, as for the three-jaw SPG in grasping state shown in Figure 1f, and the corresponding dynamic equation can be formulated as

(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32)

where $\mathbf{M}$ is the global mass matrix, $\mathbf{\Phi }=\mathbf{0}$ is the constraint equation for the SPG subjected to a fixed end constraint, and λ is the Lagrangian multiplier vector corresponding to $\mathbf{\Phi }=\mathbf{0}$. ${\mathbf{\Phi }}^{\mathrm{N}}=\mathbf{0}$ is the displacement constraint equation for the contact point between the SPG and the ball, and ${\mathit{\lambda }}^{\mathrm{N}}$ is the Lagrangian multiplier vector corresponding to ${\mathbf{\Phi }}^{\mathrm{N}}=\mathbf{0}$. ${\mathbf{Q}}_F$, ${\mathbf{Q}}_k$, ${\mathbf{Q}}_p$, ${\mathbf{Q}}_g$, ${\mathbf{Q}}_c$, ${\mathbf{Q}}_f$ severally refer to the concentrated force vector exerted on the center point of the base, the global generalized elastic force vector, pressure force vector, gravitational force vector, contact force vector and frictional force vector, which can be obtained through an assembling process from the corresponding elemental matrices or vectors and expressed as

(33) $\mathbf{M}=\sum _{\mathit{e}=1}^{n_k}{\mathbf{B}}_e^{\mathrm{T}}{\mathbf{M}}_e{\mathbf{B}}_e,{\mathbf{M}}_e={\int }_V\mathit{\rho }{\mathbf{S}}^{\mathrm{T}}\mathbf{S}\mathit{dV}$(33)

(34) $\begin{array}{rl}& {\mathbf{Q}}_k=\sum _{\mathit{e}=1}^{n_k}{\mathbf{B}}_e^{\mathrm{T}}{\mathbf{Q}}_{\mathit{ek}},{\mathbf{Q}}_p=\sum _{\mathit{e}=1}^{n_p}\mathbf{B}{{ }_e^p}^{\mathrm{T}}{\mathbf{Q}}_{\mathit{ep}},{\mathbf{Q}}_g=\sum _{\mathit{e}=1}^{n_k}{\mathbf{B}}_e^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eg}},{\mathbf{Q}}_F=\sum _{\mathit{e}=1}^{n_k}{\mathbf{B}}_e^{\mathrm{T}}{\mathbf{Q}}_{\mathit{eF}}\\ & {\mathbf{Q}}_c=\sum _{\mathit{e}=1}^{n_c}\left(\mathbf{B}{{ }_e^m}^{\mathrm{T}}{\mathbf{Q}}_{ec}^m+\mathbf{B}{{ }_e^s}^{\mathrm{T}}{\mathbf{Q}}_{ec}^s\right),{\mathbf{Q}}_f=\sum _{\mathit{e}=1}^{n_c}\left(\mathbf{B}{{ }_e^m}^{\mathrm{T}}{\mathbf{Q}}_{ef}^m+\mathbf{B}{{ }_e^s}^{\mathrm{T}}{\mathbf{Q}}_{ef}^s\right)\end{array}$(34)

where $\mathit{\rho }$ is the mass density, V is the volume of the ANCF solid element in original state, ${\mathbf{M}}_e$ is the constant symmetric element mass matrix. ${\mathbf{Q}}_{\mathit{ek}}$ is the element elastic force vector, ${\mathbf{Q}}_{\mathit{ep}}$ is the element pressure force vector, ${\mathbf{Q}}_{\mathit{eg}}$ is the element gravitational force vector, ${\mathbf{Q}}_{\mathit{eF}}$ is the element concentrated force vector, ${\mathbf{Q}}_{ec}^s,{\mathbf{Q}}_{ec}^m$ are the element contact force vectors of the slave and master surfaces, respectively, ${\mathbf{Q}}_{ef}^s,{\mathbf{Q}}_{ef}^m$ are the element tangential frictional force vectors of the slave and master surfaces, respectively. In addition, ${\mathbf{B}}_e^{\hspace{0.17em}}$, ${\mathbf{B}}_e^p$, ${\mathbf{B}}_e^m$, ${\mathbf{B}}_e^s$ are the corresponding Boolean matrices, and n_k is the number of the solid element, n_p is the number of the element subjected to air pressure, n_c is the number of the contact pairs.

Gaussian quadrature integration technique is employed to calculate the aforesaid different element force vectors and their derivatives [32, 42]. The nonlinear dynamic Equation (32)(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32) can be solved by uniting Generalized-α scheme [43] with Newton-Raphson iteration approach [44], and the time step is set as Δt = 10⁻⁴ s to guarantee convergence during the calculation process. The (k +1)th step and the (k)th step have undermentioned relation in iterative numerical calculation.

(35) $\begin{array}{rl}& \left[\begin{array}{c}{\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k+1)}\\ {\lambda }_{t+\mathrm{\Delta }t}^{(k+1)}\\ ({\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k+1)}\end{array}\right]=\left[\begin{array}{c}{\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)}\\ {\lambda }_{t+\mathrm{\Delta }t}^{(k)}\\ ({\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\end{array}\right]+\left[\begin{array}{c}\mathrm{\Delta }{\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)}\\ \mathrm{\Delta }{\lambda }_{t+\mathrm{\Delta }t}^{(k)}\\ (\mathrm{\Delta }{\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\end{array}\right]=\\ & \left[\begin{array}{c}{\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)}\\ {\lambda }_{t+\mathrm{\Delta }t}^{(k)}\\ ({\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\end{array}\right]+{\mathbf{H}}_t^{-1}\left|{ }_{\left({\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)},{\lambda }_{t+\mathrm{\Delta }t}^{(k)},({\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\right)}\right.\cdot \left[-\bar{\mathbf{\Phi }}\left({\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)},{\lambda }_{t+\mathrm{\Delta }t}^{(k)},({\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\right)\right]\end{array}$(35)

where

(36) $\begin{array}{c}{\mathbf{H}}_t=\left[\begin{array}{ccc}\frac{\mathrm{\partial }\bar{\mathbf{\Phi }}(\mathbf{q},\mathit{\lambda })}{\mathrm{\partial }\mathbf{q}}& \frac{\mathrm{\partial }\bar{\mathbf{\Phi }}(\mathbf{q},\mathit{\lambda })}{\mathrm{\partial }\mathit{\lambda }}& \frac{\mathrm{\partial }{\bar{\mathbf{\Phi }}}^{\mathrm{N}}(\mathbf{q},{\mathit{\lambda }}^{\mathrm{N}})}{\mathrm{\partial }{\mathit{\lambda }}^{\mathrm{N}}}\end{array}\right]=\\ \left[\begin{array}{ccc}\mathit{\alpha }{ }^{\mathrm{\prime }}\mathbf{M}+\mathrm{\partial }({\bar{\mathbf{\Phi }}}_{\mathbf{q}}^{\mathrm{T}}\mathit{\lambda })/\mathrm{\partial }\mathbf{q}+\mathrm{\partial }(\bar{\mathbf{\Phi }}{{ }_{\mathbf{q}}^{\mathrm{N}}}^{\mathrm{T}}{\mathit{\lambda }}^{\mathrm{N}})/\mathrm{\partial }\mathbf{q}+{\mathbf{K}}_k-{\mathbf{K}}_F-{\mathbf{K}}_p-{\mathbf{K}}_g-{\mathbf{K}}_c-{\mathbf{K}}_f& {\bar{\mathbf{\Phi }}}_{\mathbf{q}}^{\mathrm{T}}& \bar{\mathbf{\Phi }}{{ }_{\mathbf{q}}^{\mathrm{N}}}^{\mathrm{T}}\\ \overline{{\mathbf{\Phi }}_{\mathbf{q}}}& \mathbf{0}& \mathbf{0}\\ {\bar{\mathbf{\Phi }}}_{\mathbf{q}}^{\mathrm{N}}& \mathbf{0}& \mathbf{0}\end{array}\right]\end{array}$(36)

The iterative calculation process terminates when the norm of $\mathbf{\Delta }={{\left[\begin{array}{ccc}\mathrm{\Delta }{\mathbf{q}}_{t+\mathrm{\Delta }t}^{(k)}& \mathrm{\Delta }{\lambda }_{t+\mathrm{\Delta }t}^{(k)}& (\mathrm{\Delta }{\mathit{\lambda }}^{\mathrm{N}}{)}_{t+\mathrm{\Delta }t}^{(k)}\end{array}\right]}^{\mathit{\hspace{0.17em}}}}^{\mathrm{T}}$ is less than a sufficiently small positive quantity, that is, $\parallel \mathbf{\Delta }\parallel \le \mathit{\epsilon }$, where $\mathit{\epsilon }$ is the error tolerance set as $\mathit{\epsilon }$ = 10⁻⁶.

Moreover, von Mises stress distributions during grasping process of the SPG can be obtained by solving the following Equations (37)(37) ${\sigma }_s=\sqrt{\frac{{{\left({\sigma }_x-{\sigma }_y\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_y-{\sigma }_z\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_z-{\sigma }_x\right)}^{\mathit{\hspace{0.17em}}}}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}{2}}$(37) and (38(38) $\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{xy}}& {\sigma }_y& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{xz}}& {\tau }_{\mathit{yz}}& {\sigma }_z\end{array}\right]=2\left(\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}\mathit{i}\left(I_1-3\right){ }^{\mathit{i}-1}\right){\mathbf{I}}_3+\mathit{\kappa }\left(J_C-1\right)J_C{\mathbf{C}}^{-1}$(38) )

(37) ${\sigma }_s=\sqrt{\frac{{{\left({\sigma }_x-{\sigma }_y\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_y-{\sigma }_z\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_z-{\sigma }_x\right)}^{\mathit{\hspace{0.17em}}}}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}{2}}$(37)

(38) $\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{xy}}& {\sigma }_y& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{xz}}& {\tau }_{\mathit{yz}}& {\sigma }_z\end{array}\right]=2\left(\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}\mathit{i}\left(I_1-3\right){ }^{\mathit{i}-1}\right){\mathbf{I}}_3+\mathit{\kappa }\left(J_C-1\right)J_C{\mathbf{C}}^{-1}$(38)

where ${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z$ are severally first, second, and third principal stresses, and ${\tau }_{\mathit{xy}}$, ${\tau }_{\mathit{yz}}$, ${\tau }_{\mathit{zx}}$ are shear stresses.

4. Simulation analysis and experimental investigation

4.1. Simulation and experimental study of a single soft pneumatic actuator

First, we perform simulations and experimental studies on a single soft pneumatic actuator at 35 kPa pressure to verify the accuracy of the proposed model, and the theoretical results can be obtained by solving Equation (32)(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32) , (37(37) ${\sigma }_s=\sqrt{\frac{{{\left({\sigma }_x-{\sigma }_y\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_y-{\sigma }_z\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_z-{\sigma }_x\right)}^{\mathit{\hspace{0.17em}}}}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}{2}}$(37) ) and (38(38) $\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{xy}}& {\sigma }_y& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{xz}}& {\tau }_{\mathit{yz}}& {\sigma }_z\end{array}\right]=2\left(\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}\mathit{i}\left(I_1-3\right){ }^{\mathit{i}-1}\right){\mathbf{I}}_3+\mathit{\kappa }\left(J_C-1\right)J_C{\mathbf{C}}^{-1}$(38) ). As it is not easy to experimentally measure stress distribution of the soft actuator, and thus finite element simulations are carried out in ABAQUS software to validate stress distribution. Besides, we also conduct experiments on the soft actuator to further clarify the model’s accuracy. The final results are shown in Figure 5.

Figure 5. Simulations and experiments of a single soft pneumatic actuator. (a) The proposed model’s result, (b) ABAQUS result, (c) the experimental result, (d) Comparative results.

From the color bars in Figure 5 (a)–(b), it can be clearly known that the stresses essentially concentrate in the contact air chamber walls and root areas of adjacent two chambers, and the magnitude of both values is basically the same. Moreover, from the perspective of the deformed configurations, it can be seen from Figure 5 that the proposed model’s result, ABAQUS result and experimental result are basically consistent. Therefore, the model’s accuracy is validated.

4.2. Simulation analysis of the three-jaw SPG

After verifying the proposed model’s accuracy, and then the three-jaw SPG is simulated and analyzed. At the beginning, a small ball is stationary on a support platform, the three-jaw SPG is uninflated, and there is no contact interaction between the SPG and the ball. Then, the SPG is inflated uniformly until the air pressure is equal to 35 kPa at t = 1 s, and the pressure always remain unchanged, and the load curve is shown in Figure 6. The pressure should be neither too large nor too small, if it is too large, it is possible to burst the air chamber walls. On the contrary, if it is too small, the bending deformation is too small to provide a sufficient grasping force, and the force roughly reaches an appropriate state when the magnitude of the pressure is around 35 kPa. At the moment t = 2 s, the ball falls freely when the gripper grasps the ball. Then, a concentrated force (i.e. 1 N) and an initial velocity (i.e. 24 mm/s) are exerted on the entire structure at the moment t = 3 s.

Figure 6. The applied pressure load.

The detailed parameters of the SPG are listed in Table 1. The material parameters of the Yeoh model are C₁₀ = 0.2712 MPa, C₂₀ = 0.03053 MPa, C₃₀ = −0.0004013 MPa, and the incompressibility constant K = 1000 Mpa, the normal and tangential contact stiffness are set as k_n = k_τ = 4 N/mm³ [32, 39, 45]. The static and dynamic frictional coefficients between the SPG and the ball are set as μ_s = 0.4 and μ_d = 0.3, respectively. The radius and density of the ball are severally r = 20 mm and ρ = 1100 kg/m³.

Herein, neglecting viscous effect between the SPG and the ball, we only consider the contact interaction between both of them, combine a multiple-point contact method based on ANCF [32, 46] with Coulomb’s friction law [47], mesh the SPG and the ball at intervals of 1 mm, and then conduct simulation analysis on whole motion process of the soft gripper grasping the ball through numerical solution of Equation (32)(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32) .

Figure 7 shows the configurations and von Mises stress distribution of the SPG and displacement variation of the ball at different moments obtained by solving Equations (32)(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32) , (37(37) ${\sigma }_s=\sqrt{\frac{{{\left({\sigma }_x-{\sigma }_y\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_y-{\sigma }_z\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_z-{\sigma }_x\right)}^{\mathit{\hspace{0.17em}}}}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}{2}}$(37) ) and (38(38) $\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{xy}}& {\sigma }_y& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{xz}}& {\tau }_{\mathit{yz}}& {\sigma }_z\end{array}\right]=2\left(\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}\mathit{i}\left(I_1-3\right){ }^{\mathit{i}-1}\right){\mathbf{I}}_3+\mathit{\kappa }\left(J_C-1\right)J_C{\mathbf{C}}^{-1}$(38) ), which reflects the deformation, stress and displacement of the SPG in the process of grasping. At the incipient moment t = 0.0 s, the uninflated SPG hangs down naturally, and the small ball is at rest on the platform (Figure 7 (a)). As the SPG is gradually inflated, and it starts to contact with the ball at t = 0.7 s, the deformation and stress of adjacent two chambers increasingly grow large until t = 3.0 s (Figure 7 (b)–(c)), and the ball drops freely down under the action of gravitational force at the moment t = 2.0 s, whereas the overall configuration and stress distribution hardly change before t = 3.0 s. Subsequently, the pressure remains constant (i.e. p = 35 kPa) to hold the ball, and then the concentrated force 1 N is applied to the overall structure, and the SPG can lift up the ball at a constant velocity of 24 mm/s. From t = 3.0 s to t = 6.0 s, under the combined action of the air pressure and the concentrated force, the structure moves upwards by a displacement of about 70 mm (Figure 7 (c)–(f)). Since the air pressure keeps constant, the stress distribution state remains basically unchanged, and only the displacement varies.

Figure 7. The configurations and von Mises stress of the SPG and displacement variations of the ball at different moments. (a) t = 0.0 s, (b) t = 0.7 s, (c) t = 3.0 s, (d) t = 4.0 s, (e) t = 5.0 s, (f) t = 6.0 s.

One can easily observe from Figure 7 that the global stress of the undeformed entire structure is zero at initial resting state. From t = 0.7 s to t = 6.0 s, it is evident from color bars that the stresses in the contact areas and the root areas of adjacent two chambers reach the peak value, whereas the values in other areas are relatively small.

Calculation and analysis of the stress variations in these areas will facilitate the optimal design of the SPG without destroying air chambers. It is worth noting that the stresses of the base and the ball are kept at zero all the time during the grasping process, the reason is that both are completely hard bodies, they do not deform their shapes and remain unchanged. On one hand, the overall configuration of the SPG varies with time before the air pressure reaches a maximum value of 35 kPa (i.e. before t = 1.0 s). After that, it almost remains unchanged, the SPG can tightly grasp the ball and lift it to a certain height. On the other hand, the overall configuration variations and stress distribution as well as stress concentration phenomenon can offer us a distinct and straightforward view for exploring the grasping capability of the SPG, which discloses its nonlinear deformation and grasping behaviors under a certain condition.

In addition, to reveal the motion process and kinematic characteristics of grasping the ball, we carry out further study on the displacement and contact force variations of the ball. Figure 8 shows the displacement of the ball in two time periods, it can be observed that the ball is still on the support platform before t = 2.0 s, with the displacements remaining zero in three directions. Subsequently, it falls at t = 2.0 s under the action of the gravitational force. The displacement in the vertical direction (i.e. Y direction) sharply drops by about 1.3 mm in a very short time, and the ones in other two directions (i.e. X and Z directions) only fluctuate within about 0.2 mm due to frictional contact and then basically maintain constant. Starting from the moment t = 3.0 s, the displacement in Y direction varies linearly with time and reaches around 70 mm until t = 6.0 s. Meanwhile, the displacements in X and Z directions remain unchanged, and the reason lies in the fact that the applied concentrated force is equal to the sum of the gravitational forces of the SPG and the ball, which causes the entire structure to uniformly move upwards at a given initial velocity of 24 mm/s. The kinematic behaviors can be expounded in terms of displacement changes.

Figure 8. The displacement versus time of the ball in two time periods. (a) t = 0.0 s ~ 3.0 s, (b) t = 3.0 s ~ 6.0 s.

To elucidate the essence of the ball movement during grasping process, we analyze different forces exerted on the ball. Figure 9 exhibits different force variations on the ball grasped by the SPG, and it can be found that the normal contact force, tangential frictional force, and resultant force exerted on the ball along Y direction increase linearly from t = 0.7 s to t = 1.0 s when the SPG begins to be inflated and contact with the ball, and then maintain constant until t = 2.0 s. Followed by, they descend drastically at t = 2.0 s due to the ball’s free-falling motion, and instantly return to a certain value, which is induced by the fact that the gravitational force of the ball is slightly greater than the grasping force generated by the air pressure at the instant of t = 2.0 s, and the grasping force immediately exceeds self-weight of the ball after that instant. This phenomenon leads to a certain degree of relaxation of the grasping force. However, the constant air pressure rapidly restores the grasping force to one certain value and fluctuate in a small range. One can also observe from Figure 9 (a)–(b) that the normal contact force in Y direction dominates and is much greater than the tangential frictional force, and the forces in other two directions only vary in a very small scope before t = 3.0 s and remain constant after t = 3.0 s.

Figure 9. The ball subject to diverse forces. (a) Normal contact force, (b) Tangential frictional force, (c) Resultant force.

It can be known from the above comparative analysis that although the classical FEM can also characterize the bending deformation of the SPG, the proposed method can both clearly illustrate detailed configurations of the SPG and accurately describe its displacement and contact force variations at each moment from a dynamic perspective.

Defining G_SPG, G_B, F_C as the gravitational force of the SPG, the gravitational force of the ball, the concentrated force exerted on the center point of the base respectively, G_SPG and G_B can be calculated, as shown in Table 2. It can be known from Figure 9 (c) and Table 2 that in the process of grasping the ball after t = 3.0 s, the resultant force (about 0.36 N) induced by the air pressure in Y direction is nearly equal to self-weight of the ball, and thus the SPG and the ball remain relatively stationary. Moreover, after t = 3.0 s, the magnitude of the concentrated force (i.e. 1 N) approximately equals to the sum of the gravitational forces of the SPG and the ball, which satisfies the condition of uniform linear motion and causes the ball to move upward at a uniform velocity of 24 mm/s. The equilibrium relationship between different forces in Figure 9 and Table 2 enables the ball to maintain a uniform motion law, which is consistent with the displacement versus time curve in Figure 8 (b).

Table 2. Different forces in the process of uniform linear motion.

Description	Parameters	Magnitude
The gravitational force of the SPG	GSPG	0.6292 N
The gravitational force of the ball	GB	0.3685 N
The concentrated force	FC	1 N

4.3. Experimental investigation

To assess grasping capability of the SPG and examine its kinematic characteristics, an actual grasping experiment on the ball is carried out by setting up an experimental platform. Experimental results are compared with simulation ones from the theoretical modeling to draw some rules.

The SPG is fixed at the end of UR5 robotic arm (Universal Robots corporation, located in Denmark) controlled by control box. Through tuning air pressure controller (iPCU2-SMN-1120 V2.83) to regulate the pressure value in air pump (YH-04), and the deformation motion of the SPG can be adjusted. A digital force gauge (HP-20) is used to measure contact force during bending deformation of a single soft actuator, and the movement of the SPG grasping the ball can be measured by a 3D motion capture system (OptiTrack Prime 41). The entire motion process is recorded by a digital camera (Canon, EOS 1300D), and the experimental equipments are shown in Figure 10.

Figure 10. The experimental equipments.

The detailed grasping control process is as follows:

Step 1: Turn on the power supply and air pump, adjust air pressure to negative pressure, and the gripper is evacuated internally to produce reverse bending deformation;

Step 2: Control the motion of UR5 robotic arm to make the gripper gradually approach the target object;

Step 3: Adjust air pressure to positive pressure, and the gripper gradually changes from reverse bending to positive bending to adaptively contact the target object;

Step 4: Adjust air pressure to positive 35 kPa to make the gripper grasp the target object tightly;

Step 5: Control UR5 robotic arm to move upwards so as to lift the target object.

The actual total weight of the fabricated SPG and the ball is measured as 0.98 N, which is close to the theoretical value, and the error between experimental and theoretical values is probably caused by manufacturing. Then, UR5 robotic arm is programmed to rise at a uniform velocity of 24 mm/s, and the grasping states at six different moments are shown in Figure 11. The six configurations are basically an one-to-one match with that in Figure 7, which demonstrates a basic comparability between theoretical simulation results and actual experimental ones in terms of deformed configurations of the SPG.

Figure 11. Grasping the ball experiment. (a) t = 0.0 s, (b) t = 0.7 s, (c) t = 3.0 s, (d) t = 4.0 s, (e) t = 5.0 s, (f) t = 6.0 s.

For the sake of easy comparison, we can compare the experiment displacements with simulation results of the ball in Y direction. It can be known from Figure 12 that the simulation results coincide with the experimental ones, which is consistent with the kinematic law of the ball before and after being grasped.

Figure 12. The displacement of the ball in Y direction in two time periods. (a) t = 0.0 s ~ 3.0 s, (b) t = 3.0 s ~ 6.0 s.

Additionally, we need to consider another situation where there is no frictional contact between the SPG and the ball, and there is also no self-contact for the SPG itself. Similarly, we can conduct simulation and obtain the deformation configuration and von Mises stress variations of the SPG at t = 2.5 s in the case of no frictional contact and self-contact, as shown in the overall view and partial enlarged cross-sectional view in Figure 13. It can be obviously observed from the white oval area in Figure 13 that an interpenetration phenomenon occurs between adjacent air chambers, and between the SPG and the ball, which indicates that there is no contact interaction between them. They simply undergo their own configuration and displacement variations, which is inconsistent with the practical grasping contact situation.

Figure 13. The deformation configuration and von Mises stress distribution of the SPG at t = 2.5 s.

Furthermore, we can obtain the motion changes of the ball and von Mises stress variations of the SPG by working out Equations (32)(32) $\overline{\Phi }=\left[\begin{array}{l}\mathit{M}\stackrel{\mathrm{..}}{q}+{\Phi }_q^{\text{T}}\mathit{\lambda }+{\Phi }_q^{\text{NT}}{\lambda }^{\text{N}}+Q_k-Q_F-Q_p-Q_g-Q_c-Q_f\\ \mathit{\Phi }\\ {\Phi }^N\end{array}\right]=0$(32) , (37(37) ${\sigma }_s=\sqrt{\frac{{{\left({\sigma }_x-{\sigma }_y\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_y-{\sigma }_z\right)}^{\mathit{\hspace{0.17em}}}}^2+{{\left({\sigma }_z-{\sigma }_x\right)}^{\mathit{\hspace{0.17em}}}}^2+6\left({\tau }_{xy}^2+{\tau }_{yz}^2+{\tau }_{zx}^2\right)}{2}}$(37) ) and (38(38) $\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{\mathit{xy}}& {\tau }_{\mathit{xz}}\\ {\tau }_{\mathit{xy}}& {\sigma }_y& {\tau }_{\mathit{yz}}\\ {\tau }_{\mathit{xz}}& {\tau }_{\mathit{yz}}& {\sigma }_z\end{array}\right]=2\left(\sum _{\mathit{i}=1}^3{C}_{\mathit{i}0}\mathit{i}\left(I_1-3\right){ }^{\mathit{i}-1}\right){\mathbf{I}}_3+\mathit{\kappa }\left(J_C-1\right)J_C{\mathbf{C}}^{-1}$(38) ), as shown in Figure 14. Initially, the ball is stationary (Figure 14 (a)), starting from t = 2.0 s, the ball is in a free-falling state (Figure 14 (b)-(c)), and its displacement versus time curve is half that of the parabola (Figure 14 (d)), which conforms to the law of free-falling motion. However, as is well known, it is almost impossible to conduct grasping operation without frictional contact, as it obviously goes against the motion law of the actual grasping objects.

Figure 14. The motion states of the ball and von Mises stress distribution of the SPG at (a) t = 0.0 s, (b) t = 2.0 s, (c) t = 3.0 s and (d) the displacements of the ball.

In summary, the dynamic modeling method provides a workable solution for accurately analyzing nonlinear deformation and configuration variations of the SPG and simulating contact interaction to avoid interpenetration problem when grasping target objects, which demonstrates the crucial role and influence of contact interaction on grasping manipulation of target objects.

Although the proposed method has lower computational efficiency compared to finite element software, it can stably and accurately predict the output results under external interference when considering contact interaction, including configuration variations, von Mises stress distribution and contact forces, which is instructive and enlightening for actual grasping.

5. Conclusions and further work

In this study, within the framework of ANCF, in consideration of contact interaction between soft gripper and target objects, we present a nonlinear dynamic modeling approach for a three-jaw SPG, whose grasping capability and motion behaviors can be revealed and predicted effectually.

The key contributions of this study are summarized as follows: (i) Creatively combining Coulomb’s friction law with a multiple-point contact method to establish a nonlinear dynamic model for illuminating the SPG’ grasping characteristics, several nonlinear forces and their derivatives are deduced, and the established nonlinear dynamic equations are solved using Generalized-α scheme and Newton-Raphson iterative algorithm; (ii) The proposed nonlinear dynamic approach can not only accurately describe motion displacement, but also reveal the SPG’s grasping configuration, stress distribution and contact forces variations; (iii) Revealing the pivotal role and influence of contact interaction in circumventing interpenetration problem and grasping operation.

In the future, we will be engaged in developing some new computing strategies with the help of data-driven method to improve computational efficiency and balance calculation accuracy. Additionally, efforts will be committed to implanting soft force sensors into SPGs to measure grasping force changes, and exploring the effect of different air pressures on nonlinear deformation behaviors to provide real-time feedback and accurate control, which may inspire potential applications of soft grippers in sophisticated scenes such as deepwater operations and surgical treatments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was supported by Natural Science Foundation of Zhejiang Province (Grant No. LQ22A020003) and National Natural Science Foundation of China (Grant No. 52075499) for which all authors are grateful. This research was also supported by China Scholarship Council.