Mechanical Deformation Analysis of a Flexible Finger in Terms of an Improved ANCF Plate Element

Abstract

In recent years, flexible continuum robots have been substantially developed. Absolute nodal coordinates formulation (ANCF) gives a feasible path for simulating the behavior of flexible robots. However, the model of finger-shaped robots is often regarded as a cylinder and characterized by a beam element. Obviously, this is short of characterizing the geometrical feature of fingers in detail, especially under bending conditions. Additionally, for the lower-order plate element, it is hard to characterize the bending behavior of the flexible finger due to fewer nodes; a higher-order plate element often requires an extremely long computing time. In this work, an improved ANCF lower-order plate element is used to increase the accuracy of the Yeoh model and characterize the geometrical structure of silicone rubber fingers, taking into particular consideration the effect of volume locks and multi-body system constraints. Since it is a kind of lower-order plate element, essentially, the computing time is nearly the same as that of conventional lower-order plate elements. The validity of this model was verified by comparing it with the results of the published reference. The flexible finger, manufactured using silicone rubber, is characterized by the novel ANCF lower-order plate element, whereby its mechanical deformation and bending behavior are simulated both efficiently and accurately. Compared to the ANCF beam element, conventional lower-order plate element, and higher-order plate element, the novel plate element in this paper characterizes the external contour of the finger better, reflects bending behavior more realistically, and converges in less computing time.

1. Introduction

As an indispensable part of mechanical grippers and robots, the structure of silicone rubber flexible fingers has taken many forms and changed in complexity under the effort of many scholars. Its origin can be traced back to previously pneumatic artificial muscles, which were designed to simulate human muscles and consisted of silicone rubber cylinders and outer woven nets. Since then, the property of silicone rubber has improved, and the volume of artificial muscles has shrunk gradually and has moved towards soft biological tissue because soft biological tissue has unique advantages compared to hard structures, such as strong environmental adaptability, flexibility, and being harmless to the human body. The silicone rubber flexible finger has such characteristics and uses them for the grasping process; its excellent packing skill is superior to traditional rigid mechanical grippers in food handling, medical service robots, and so on.

Silicon rubber materials are classified as hyper-elastic materials in terms of mechanical properties, the attribution of which should be characterized by hyper-elastic constitutive equations. The Neo-Hookean model is one of the most commonly used constitutive models of hyper-elastic materials; however, it is a simple model that expresses stain energy, and it only has good agreement under small strains; the computation will be less accurate as the strain energy increases. Mooney researched the deformation of hyper-elastic materials theoretically and experimentally and then raised the stain energy model of the first and second invariants for the right Cauchy–Green deformation tensor [1]. The Yeoh model is a branch of Mooney’s model; its advantage is being more suitable for calculating larger stress and strain problems as it can characterize the deformation behavior of silicone rubber material more accurately [2].

Absolute nodal coordinates formulation (ANCF) is an emerging finite element method (FEM). Its node coordinates include slope coordinates, the element mass matrix is constant, and the stiffness matrix and nodal coordinates are highly correlated. The theory of the ANCF method does not contain the hypothesis of small turning angle motions [3,4,5,6,7], which is the most significant advantage compared to the classic FEM; hence, it is suitable for calculating large deformation problems such as modern flexible fingers or other flexible continuum robots [8,9,10,11,12,13,14].

Jung discussed the accuracy of three hyper-elastic constitutive models under ANCF beam elements in modeling silicone rubber, and the result shows that the Yeoh model is more accurate [15]. Pappalardo et al. proposed a new fully parametrized plate finite element based on ANCF kinematic description that fills the theoretical gap of the large rotation vector formulations [16]. Melly et al. obtained the residual strain energy density from equibiaxial loading and fitted it to a term with dependence on the second invariant; then, they added the term to the original expression, whereby modifying Yeoh’s model for improving the prediction of vulcanized rubber and thermoplastic elastomers [17]. A rational absolute nodal coordinate formulation (RANCF) thin plate element was developed by Pappalardo; it is used in the accurate geometric modeling and analysis of flexible continuum bodies with complex geometrical shapes [18]. The flexible body may have an inaccurate rigid characteristic when the incompressible material model encounters bending deformation. Namely, the model of the hyper-elastic incompressible material is unable to converge to a correct result, the phenomenon of which is named volume lock. In order to solve this problem, Bayat et al. presented a locking-free element formulation based on reduced integration; physically-based hourglass stabilization (Q1SP) was coupled for the first time with the DG framework, leading to a DG variant with very good convergence properties [19]. Orzechowski et al. established an efficient method of modeling nonlinear nearly incompressible materials with polynomial Mooney–Rivlin models and a volumetric energy penalty function in the ANCF framework [20]. This method was dedicated to the examination of several ANCF fully parameterized beam elements under an incompressible regime. Xu et al. proposed a higher-order plate element formulation with quadratic interpolation in the transverse direction, which can not only alleviate volumetric locking but also improve accuracy in the simulation of large bending deformations compared to improved lower-order plate elements with the selectively reduced integration method [21].

In the published references, models of finger-shaped robots are often regarded as a cylinder and characterized by a beam element. Obviously, this is short of characterizing the geometrical feature of the fingers in detail, especially under bending conditions. Additionally, for lower-order plate elements, it is hard to characterize the bending behavior of the flexible finger due to fewer nodes; higher-order plate elements often require extremely long computing times.

In this work, an improved ANCF lower-order plate element is used to increase the accuracy of the Yeoh model and characterize the geometrical structure of the silicone rubber finger, taking into particular consideration the effect of volume locks and multi-body system constraints. Since it is a kind of lower-order plate element, essentially, the computing time is nearly the same as that of conventional lower-order plate elements. The validity of this model was verified by comparing it with the results of the published reference. The flexible finger, manufactured using silicone rubber, is characterized by the novel ANCF lower-order plate element, whereby its mechanical deformation and bending behavior are simulated both efficiently and accurately.

This paper is organized as follows. After the introduction, an improved Yeoh model based on the ANCF plate element is introduced in Section 2. In Section 3, the relevant examples’ verification is conducted and analyzed. On the basis of the above-mentioned theory, the computation of a silicone rubber finger is presented in Section 4. In the last section, we give a short conclusion.

2. Theoretical Foundation and Modeling

2.1. ANCF Plate Element Theory

2.1.1. Kinematics Characterization

A three-dimensional ANCF plate element with four nodes is shown in Figure 1, in which a, b, and t represent the length, width, and depth, respectively. Each node contains 12 node coordinates and, at any point, includes 12 degrees of freedom (DOFs), which are written in the form of a vector $e_i$,

$$e_i={\left[\begin{array}{cccc}{r}_i^{\mathrm{T}}& \frac{\partial r_i^{\mathrm{T}}}{\partial x}& \frac{\partial r_i^{\mathrm{T}}}{\partial y}& \frac{\partial r_i^{\mathrm{T}}}{\partial z}\end{array}\right]}^{\mathrm{T}}$$

(1)

where $r_i$ is the global coordinate vector of the node (A, B, C, or D); $\partial r_i/\partial \alpha$ ($\alpha =x,y,z$) is the partial derivative of the global coordinate vector to the element coordinate.

The coordinates of a single element are given in the form of a vector:

$$e={\left[\begin{array}{cccc}{e}_A^{\mathrm{T}}& e_B^{\mathrm{T}}& e_C^{\mathrm{T}}& e_D^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$$

(2)

The coordinates of any point in the element coordinate system are $x$, $y$, $z$, and its position coordinate $r={\left[\mathrm{X},\mathrm{Y},\mathrm{Z}\right]}^{\mathrm{T}}$ in the global coordinate system is given by the following position interpolation function:

$$r=S(x,y,z)e$$

(3)

where $r$ is the global coordinate vector of any point in the plate element, $S$ is the shape function matrix, and $e$ is the coordinate vector of the element, which is given by Formula (2).

$S$ is only related to the element coordinate. Different shape function matrices and node coordinate vectors are multiplied to obtain different global coordinate interpolation functions; the shape function of a typical three-dimension four-node plate element is as follows,

$$\begin{array}{l}S=[\begin{array}{cccc}{S}_{1}I& S_{2}I& S_{3}I& \begin{array}{cccc}{S}_{4}I& S_{5}I& S_{6}I& \begin{array}{cc}{S}_{7}I& S_{8}I\end{array}\end{array}\end{array}\\ \begin{array}{cccc}{S}_{9}I& S_{10}I& S_{11}I& \begin{array}{cccc}{S}_{12}I& S_{13}I& S_{14}I& \begin{array}{cc}{S}_{15}I& S_{16}I]\end{array}\end{array}\end{array}\end{array}$$

(4)

2.1.2. The Mass Matrix, Elastic Force, and External Force of an Element

The cell mass matrix can be obtained by the kinetic energy formula

$$T=\frac{1}{2}{\displaystyle \underset{V}{\int }{\rho }_0}{\dot{r}}^{\mathrm{T}}\dot{r}dV$$

(5)

where ${\rho }_0$ is the density at the beginning of the structure; V is the volume at the initial moment. We substitute Formula (3) into the above formula to obtain the element mass matrix.

$$\mathrm{M}={{\displaystyle \underset{V}{\int }{\rho }_{0}S}}^{\mathrm{T}}SdV$$

(6)

Because the shape function of the element is not related to time, in the beginning, the density and size of the element are constant, so the result of Formula (6) is a constant matrix that does not change by time. The node coordinate of the plate element is a $48\times 1$ column vector, the shape function is a $3\times 48$ matrix, and the dimension of the mass matrix of the element is $48\times 48$.

The elastic force can be obtained from stain energy, which can be characterized as:

$$U=\frac{1}{2}{\displaystyle \underset{V}{\int }{\mathit{\epsilon }}^{\mathrm{T}}}E\mathit{\epsilon }dV$$

(7)

where $\mathit{\epsilon }$ is the vector form of the strain tensor, and $E$ is the elastic coefficient matrix.

The strain tensor ${\mathit{\epsilon }}_m$ is obtained by using the right Cauchy–Green deformation tensor [22],

$${\mathit{\epsilon }}_m=\frac{1}{2}(J^{\mathrm{T}}J-I)$$

(8)

$J$ is the displacement gradient matrix.

$$J=\frac{\partial r}{\partial r_0}$$

(9)

$r$ and $r_0$ represent the coordinate vector of the element at this moment and the coordinate matrix of the element at the initial moment, respectively.

Formula (7) means that the elastic force Q_s is computed by the partial differential of the nodal coordinate vector due to stain energy.

$$Q_s={\left(\frac{\partial U}{\partial e}\right)}^{\mathrm{T}}$$

(10)

Any point on the element is subjected to a concentrated force $F$, and the virtual work is done by this force,

$$\delta W_{\mathrm{e}}=F^{\mathrm{T}}\delta r=F^{\mathrm{T}}S\delta e=Q_e^{\mathrm{T}}\delta e$$

(11)

where $r$ is the global coordinate vector, which needs to be substituted into the coordinate values x, y, z of the element. $Q_e$ is the generalized external force corresponding to the concentrated force of the element

$$Q_e=S^{\mathrm{T}}F$$

(12)

2.1.3. Dynamics Equations

The previous parts are the theory of a single board element according to the virtual work; the dynamic equation of the entire system that can be gained under the ANCF undamped and unconstrained dynamic equation is

$$M{\ddot{e}}_s+Q_s-Q_e=0$$

(13)

where $M$ is assembled from the mass matrix of the element and is not changed with time; $e_s$ is the system coordinates, which are assembled from element coordinates; $Q_s$ is the elastic force vector of the system, which is assembled from the elastic force of the element. The elastic force and the node coordinates show extreme nonlinearity with time change; $Q_e$ is the generalized external force vector of the system [23].

For a system containing constraints and damping, the dynamic equation is

$$M{\ddot{e}}_s+C{\dot{e}}_s+{\mathsf{\Phi }}_{e_s}^{\mathrm{T}}\lambda +Q_s-Q_e=0$$

(14)

where C is the system damping matrix; $\mathsf{\Phi }$ is the constraint equation; $\lambda$ is the Lagrangian multiplier vector.

2.2. An Improved Hyper-Elastic Constitutive Model Based on ANCF

2.2.1. A Constitutive Model Based on ANCF

In the theory of continuum mechanics, the three invariants of the right Cauchy–Green deformation tensor $C_r=J^{\mathrm{T}}J$ are, respectively,

$$I_1=\mathrm{tr}(C_r)$$

(15)

$$I_2=\frac{1}{2}\left[\mathrm{tr}{(C_r)}^2-\mathrm{tr}(C_r^2)\right]$$

(16)

$$I_3=\det (C_r)$$

(17)

In the above three formulas, $I_1$, $I_2$ and $I_3$ are, respectively, the first invariant, the second invariant, and the third invariant of the right Cauchy–Green deformation tensor. When the material is incompressible, $J=\det (J)=1$, or equivalently, $I_3=J^2=1$. Therefore, the strain energy of incompressible materials depends only on $I_1$ and $I_2$.

The Yeoh model is a function of the first invariant of the right Cauchy–Green deformation tensor, and its model form is:

$$U_y={\displaystyle \sum _{i=1}^3{\mu }_{i0}{(I_1-3)}^i}$$

(18)

In the above formula, ${\mu }_{i0}$ is the material constant evaluated by experiments. Although the first invariant was used in the Yeoh model, its order reaches the third order. At the same time, there are three material constants that characterize material properties that have been included in this model, and the content of expression is richer.

The model expressed by Equation (18) assumes that the material is completely incompressible, which satisfies $J=\det (J)=1$. This condition needs to be added to the equation of the system; the usual methods include the Lagrangian multiplier method and the penalty method. The Lagrangian multiplier method is a constraint method that multiplies the condition of volume incompressibility by a Lagrangian multiplier and substitutes it into the system equation to achieve volume control. Compared with the Lagrangian multiplier method, the penalty method is simpler and more effective here, so the penalty method was chosen for this article to achieve volume constraints. Under the penalty method, additional strain energy density will be added to the constitutive model,

$$U_p=\frac{1}{2}k{(J-1)}^2$$

(19)

where $U_p$—the volume penalty equation; $k$—the volume penalty number, which can represent the bulk modulus of the materials. A reasonable value of $k$ is usually selected to ensure the incompressibility of the volume and set to 1000 MPa. The reason is that an excessively large value of $k$ may lead to numerical problems during the process of computation. It should be pointed out that the volumetric energy penalty equation was used to make the materials nearly incomprehensible rather than absolutely incomprehensible.

The full formula of the Yeoh model, using the penalty equation, is

$${\overline{U}}_y=U_y+U_p={\displaystyle \sum _{i=1}^3{\mu }_{i0}{(I_1-3)}^i+\frac{1}{2}}k{(J-1)}^2$$

(20)

According to the conclusion obtained by Orzechowski [20], volume behavior was introduced to the above equation, which may result in a numerical problem during the process of computation, so it is necessary to separate the volume and partial derivative components. In order to decouple, the invariant of the right Cauchy–Green deformation tensor was not directly used in this constitutive equation, which was replaced by the following form:

$${\overline{I}}_1=J^{-2/3}{I}_1$$

(21)

The above variation was used in this article, which avoided the numerical problem during the process of finding a solution. The new constitutive model after the replacement was as follows:

$${\overline{U}}_y=U_y+U_p={\displaystyle \sum _{i=1}^3{\mu }_{i0}{({\overline{I}}_1-3)}^i+\frac{1}{2}}k{(J-1)}^2$$

(22)

The above formula is a general strain energy density formula; the volume integration of each element is needed to obtain the strain energy of the element. Therefore, with the improved Yeoh model, the strain energy of the element can be expressed as

$${\overline{U}}_Y={\displaystyle {\int }_V\left(U_y+U_p\right)dV}={\displaystyle {\int }_V\left({\displaystyle \sum _{i=1}^3{\mu }_{i0}{({\overline{I}}_1-3)}^i+\frac{1}{2}}k{(J-1)}^2\right)dV}$$

(23)

Hence, the immutable volume of the Yeoh model was established through the penalty equation, and it was further improved by separating the volume and partial derivation components of deformation.

According to Formula (10), the elastic force of the improved Yeoh model with ANCF is shown as follows:

$${\overline{Q}}_{KY}={\left(\frac{\partial {\overline{U}}_Y}{\partial e}\right)}^{\mathrm{T}}={\displaystyle {\int }_V{\left(\frac{\partial U_y}{\partial e}+\frac{\partial U_p}{\partial e}\right)}^{\mathrm{T}}}dV$$

(24)

2.2.2. The Elimination of Volume Lock

Many researchers have mentioned that ANCF is affected by the lock phenomenon of the element. The lock phenomenon of the continuum mechanics method is particularly obvious, based on which the wrong rigidity characteristics are shown in the bending process of flexible bodies. The influence of this phenomenon is even greater when it comes to the incompressible material model. The super-elastic incompressible material model cannot converge to the correct result. This problem is called the volume lock phenomenon in the reference.

Numerical integration is usually used in the above elastic force formula derived from the Yeoh model. In finite element analysis, the Gaussian integral is the most commonly used numerical integration algorithm. In the ANCF elastic material model and many other finite element methods, selective reduction integration (short for SRI) is often used to prevent the occurrence of volume lock. The integration of strain energy density is divided into two parts with the point of the above method. The first part does not consider the influence of volume locks and performs a complete integration. The second part takes the impact of the rate of volume locks into consideration when using SRI.

From Formula (23), the equation of strain energy under the Yeoh model can simply be split. By removing the volumetric energy penalty equation that imposes constraints on the volume, the first part does not consider the behavior of the volume; hence, it can be completely integrated. The influence of volume is introduced to the second part through the volume penalty equation; under the plate element, it is assumed that there is no thickness, and the thickness is directly multiplied into the equation for the plane integral. The computation of Formula (23) of the strain energy equation is calculated as follows:

$${\overline{U}}_Y^{SRI}={\displaystyle {\int }_V{U}_y}dV+t{\displaystyle {\int }_S{U}_{p}dS}$$

(25)

In the above two formulas, SRI represents selective reduction integration, and t and S represent the thickness and area of the plate element, respectively. The above formula shows that the volumetric energy penalty equation is only evaluated on the neutral plane of the plate element.

2.3. The Modeling of Length Constraint

2.3.1. The Constraints of a Multi-Body System

The above system dynamic equations are unconstrained equations. There will be some constraints in the actual calculation process. The Lagrange multiplier method is usually used to compute the constraint equations of system equations. At this time, the system dynamic equations become the following form,

$$\{\begin{array}{c}M{\ddot{e}}_s+{\mathsf{\Phi }}_{e_s}^{\mathrm{T}}\lambda +Q_s-Q_e=0\\ \mathsf{\Phi }\left(e_s,t\right)=0\end{array}$$

(26)

where $\mathsf{\Phi }$ is the constraints equation, ${\mathsf{\Phi }}_{e_s}$ is the partial derivative matrix of the constraints equations with respect to the system coordinates, and $\lambda$ is the Lagrange multiplier.

Take the second derivation of the constraint equation system with respect to time t and in conjunction with the first formula in the above formula.

$$\left[\begin{array}{cc}M& {\mathsf{\Phi }}_{e_s}\\ {\mathsf{\Phi }}_{e_s}& 0\end{array}\right]\left[\begin{array}{c}{\ddot{e}}_s\\ \lambda \end{array}\right]=\left[\begin{array}{c}{Q}_e-Q_s\\ Q_d\end{array}\right]$$

(27)

where $Q_d$ is the result of the second partial derivation of the constraint equation $\mathsf{\Phi }$ with respect to time, which has the following form:

$$Q_d={\mathsf{\Phi }}_{e_s}{\ddot{e}}_s=-{\mathsf{\Phi }}_{tt}-2{\mathsf{\Phi }}_{e_{s}t}{\dot{e}}_s-{\left({\mathsf{\Phi }}_{e_s}{\dot{e}}_s\right)}_{e_s}{\dot{e}}_s$$

(28)

Hence, the conventional ordinary differential solution to Equation (27) can calculate the system response with constraints. The explicit solution process will not be discussed in detail here.

2.3.2. Length Constraint

Certain dimensionalities need to be constrained in order to be consistent with the actual condition, such as the skeleton of a flexible finger. In this subsection, a length constraint method is proposed to be applied to the plate element and proven to be both reliable and feasible in the section of theoretical verification.

As plotted in Figure 2, the length, width, and thickness of the initially undeformed single plate element are a, b, and t. Taking the AB edge as an example, if the length of the edge is unchanged during the deformation process, $l_{AB}$ should always be equal to $a$. In the process of motion deformation, the position vector of any point on the $AB$ edge in the global coordinate system can be obtained according to Formula (3):

$$r=S(x,y=0,z=0)e$$

(29)

According to the above formula, at the current moment, X, Y, and Z coordinates of any point on the AB edge in the global coordinate system are related to the element coordinate x, namely, $X=X(x)$, $Y=Y(x)$, $Z=Z(x)$. The AB curve calculated by interpolation is smooth and has continuous derivation due to ANCF. According to geometric theory, if a smooth or piecewise smooth curve has a continuous derivation, the length of the curve can be obtained by integrating the curve of the length of the arc.

$$l_{AB}={\displaystyle {\int }_{L}ds={\displaystyle {\int }_0^a\sqrt{{\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2}dx}}$$

(30)

Therefore, the constraint equation at this time can be written as:

$$\mathsf{\Phi }(e,t)=l_{AB}-a={\displaystyle {\int }_0^a\sqrt{{\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2}dx}-a=0$$

(31)

The above equation is substituted into the system of dynamics equation through Lagrangian multipliers as a constraint equation, which can make the AB edge inextensible. However, in the actual calculation process, the content of the root number in Formula (31) is too complicated to open the root number, so it is impossible to calculate the AB length at the current moment. In order to solve this problem, it is unavoidable to seek the AB length.

According to Formula (31), the coordinate of the material point of the AB segment has been determined by the parameter equation: $X=X(x)Y=Y(x)Z=Z(x)$. Hence, the length of AB is the value of the integral of the square root of the formula ${\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2$. When the length of AB is a fixed value, the integration of ${\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2$ will also be fixed. Taking ${\displaystyle {\int }_0^a{\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2}dx$ as a fixed value has the same effect as setting ${\displaystyle {\int }_0^a{\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^2}dx$ as a fixed value. Additionally, the former has the capacity to integrate.

The constraint equation can be calculated directly by removing the root sign, whereby the length of the plate element is unchanged.

$$\mathsf{\Phi }(e,t)=l_{AB}-a={\displaystyle {\int }_0^a{\left[X^{\prime }(x)\right]}^2+{\left[Y^{\prime }(x)\right]}^2+{\left[Z^{\prime }(x)\right]}^{2}dx}-a=0$$

(32)

3. Theoretical Verification

In order to verify the accuracy of the hyper-elastic constitutive model and the validity of volume lock and length constraint, some comparing calculations were conducted in this section.

3.1. Verification of Volume Lock

The simulation example is a silicone rubber plate; the structure is shown in Figure 2: one side of the plate is fixed and the other side is set free. The relevant computing parameters are listed in

Table 1. Geometrical and material parameters of a silicone rubber plate in this example.

Name	Symbol	Value	Unit
Length	a	0.18	m
Width	b	0.16	m
Thickness	t	0.04	m
Damping coefficient	c	1.5	N·s/m
Bulk modulus	k	1000	MPa
Material constant	${\mu }_{10}$	0.2712	MPa
Material constant	${\mu }_{20}$	0.03053	MPa
Material constant	${\mu }_{30}$	−0.0004013	MPa

. A dynamics simulation is carried out to calculate the movement process of the flat pendulum under the action of gravity (along the negative Z-axis), and the total calculation time is 0.4 s. The dynamics equation is established according to Equation (14), in which the damping matrix is distributed according to Rayleigh damping. The detailed explanation of each parameter is referred to in Xu’s work [21]. It is noted that ${\mu }_{10}$, ${\mu }_{20}$, and ${\mu }_{30}$ are material constants used to establish strain energy function, which requires ${\mu }_{20}$ or ${\mu }_{30}$ to be negative to meet the condition that the strain energy density is equal to or greater than zero and increases with deformation monotonically [24]. These three parameters can be obtained by tensile test [15].

In order to verify whether the SRI method eliminates the volume lock effectively, Formula (23) is used to calculate the model that does not include SRI in Figure 3; SRI is calculated using Formula (25). In the model, P is the midpoint of its side; the changes in X and Z coordinates of point P over time are shown in Figure 4 due to gravity.

It can be seen from Figure 4 that for the unreduced integration model, point P on the silicone rubber plate has no significant displacement in the X coordinate, and the maximum displacement on the Z-axis does not have the capacity to reduce half of the integration model, which moves to its limit in the vicinity of 0.15 s and then rebounds upward. Figure 5 shows the Yeoh model’s non-selectively reduced internal plate at t = 0.2 s. It can be seen that under the influence of the volume lock, the edge of the plate is tilted upward, and the local edge is even higher than the XY horizontal plane; thereby, the deformation is incorrect. After adopting SRI, the displacement of point P is significantly increased in both the X-axis and Z-axis directions, and the results of motion and unreduced integration are completely different. From the above comparison, it can be concluded that under the influence of volume locks, the silicone rubber plate cannot be swung freely, and the SRI method can effectively eliminate volume locks and avoid obtaining incorrect results.

3.2. Verification of an Improved Yeoh Model

On the basis of eliminating the effect of volume lock, the accuracy of the constitutive model is verified further. A comparison is conducted between the results of the improved Yeoh model and the experimental results in the published reference [21]. It needs to be pointed out that the comparison here is the result of improving the low-order plate element (conventional plate element by using the SRI method) rather than the results of the high-order plate element. The simulating parameters are the same as those in Section 3.1.

The comparison among the simulation results of point P in both X and Z coordinates, the experimental results, and conventional Yeoh model results is shown in Figure 6 and Figure 7. In Figure 6, the two models in the first 0.25 s have good agreement with the experimental results. In the interval from 0.25 s to 0.4 s, there is a deviation between the two models, and the result obtained by the improved Yeoh model is slightly worse.

In Figure 7, in the first 0.18 s, the results of the two models are close to the experimental results. To be specific, the improved Yeoh model seems better than the initial model. In the interval from 0.18 s to 0.3 s, the results of the improved model coincide with those obtained by the experiments. Their motion’s lowest points coincide with each other and are smaller than that of the initial Yeoh model. In the last 0.1 s, the results in terms of the two models fail to match the experimental results. What is more, the improved model shows a larger fluctuation than the others at around 0.35 s. On the whole, the improved Yeoh model that separates the volume components is in good agreement with the experiment.

Figure 8 shows the shape of the silicone rubber plate under the improved Yeoh model at various times. It can be seen from Figure 6 and Figure 7 that the plate is still swinging freely at 0.1 s, and it drops to the lowest point at 0.2 s, which is the smallest value in the Z coordinate in Figure 7. At this moment, a certain amount of elastic and kinetic energy has been accumulated on the plate, and it then continues to swing backward. In the 0.2 to 0.3 s movement, the improved Yeoh model moves faster and then slows down smoothly after 0.3 s. It can be seen from the whole process that the law of motion and shape is in accord with the theoretical expectations.

According to the whole process of the free swing of the silicone rubber plate, the Yeoh model, without any improvement, has a large deviation from the actual results. The method proposed by the published reference [20] is used to remove the coupling behavior of the volume component and the partial conductance component during the deformation process and replace the first invariant in the constitutive equation with ${\overline{I}}_1=J^{-2/3}{I}_1$. The result calculated by this improved method is very close to the experimental result, which is better than the Yeoh model. Meanwhile, computing accuracy can also be guaranteed.

3.3. Verification of Length Constraint

In this subsection, the verification calculation example is designed to illustrate the feasibility of the length constraint. As plotted in Figure 9, there is a hollow thin-walled cylinder with the specific dimensions shown in

Table 2. Parameters of the hollow cylinder.

Name	Symbol	Value	Unit
Radius	r	0.2	m
Length	L	0.4	m
Thickness	t	0.02	m
Density	$\rho$	7200	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
Elastic modulus	E	1 × 106	MPa
Poisson’s ratio	$\mu$	0.3	-

. In order to observe the deformation of the thin-walled cylinder with great flexibility, both ends of the cylinder are fixed. It should be noted that the purpose of this subsection is to verify the feasibility of the length constraint; therefore, the parameters in this case (Table 2) are totally different from those in Table 1.

Air pressure is applied to the interior to make it expand. A constraint condition of constant length is applied to the middle line where point A is located; thus, it only suffers the air load; without expansion, the left and right parts expand and deform separately. The undamped dynamic equation is adopted to establish the dynamic equation according to Formula (26), and the pressure load is applied under ANCF. The applied air pressure increases linearly with time: $\mathrm{P}=100\cdot \mathrm{t}$ kPa; the loading time is 0.5 s, which means the final value of the air pressure is 50 kPa. The deformation effect of a hollow cylinder at different moments is shown in Figure 10.

It can be seen from Figure 10 that the expanded deformation of the hollow thin-walled cylinder structure increases with air pressure. The air pressure is slight at first, and the deformation is not obvious. Significant deformation can be found at 30 kPa; due to the constraint, there is no bulge in the middle section, and symmetrical expansion and deformation only occur on the left and right sides of the constrained part; this deformation is in line with the expected deformation effect. In order to illustrate the effectiveness of the constraints in detail, the change curves in the Z-axis coordinates of points A and B in Figure 9 over time are shown in Figure 11.

In Figure 11, it is clear that the Z-axis coordinates of points A and B both become larger. At 0.5 s, when 50 kPa air pressure is loaded, the expanded deformation of point A and point B is 0.0025 and 0.035, respectively. The expanded deformation of point A does not mean that the constraint is not effective; it can be explained that the system dynamics equation is solved by a numerical algorithm, which will cause errors in the calculation process, and the cylinder itself has an increasing tendency to deform under the driving force. The errors, with an increasing trend, are accumulated during the iterative process, so the obtained results augment gradually. Comparing the small expansion of point A with the displacement of point B, the displacement of the two points is not in the same order of magnitude, which implies that the constraint effect is very good. It is feasible to use Equation (32) instead of (31) as the constraint equation, and the inextensible constraint is realized using this method. It should be pointed out that the result of Figure 11 shows a wave-like rising trend; this is because the system adopts an undamped dynamic equation without energy dissipation, and each moment is not precisely static; the external forces that augment stepwise cannot balance the elastic force of the system immediately, and it is reasonable to expect a resulting antagonistic rise.

4. The Calculation of a Silicone Rubber Finger

4.1. The Structure of a Silicone Rubber Finger

In this section, a silicone rubber flexible finger is designed, and its structure is shown in Figure 12. The flexible finger is divided into four sections. The bottom and middle sections are the same flexible cavities made of silicone rubber material. Between the two flexible cavities is a rigid bone ring supported by a more rigid material. The same silicon rubber material is used for physical packaging at the top. Except for the top, the entire flexible finger is hollow, and the thickness of the surface is uniform. The pneumatic flexible finger can be given different sizes, according to actual conditions, to meet the needs of different purposes. Some additional constraints can be added to improve its performance in the future. The deformation of the flexible finger under air pressure will be computed in this section. It should be noted that in order to improve the overall simulation rate, the overall size of the flexible finger is deliberately enlarged. The relevant deformation is still proportional to its size, so the larger size result can be used to derive the actual deformation result. The dimensions of a, b, r, and t in Figure 12 are listed in

Table 3. Parameters in Figure 12.

Name	Symbol	Value	Unit
Radius	r	0.1	m
Thickness	t	0.02	m
Length	a	0.2	m
Width	b	0.05	m

, and the relevant material parameters are the same as in Table 1.

4.2. Single Cavity Deformation Analysis

As mentioned above, this pneumatic flexible finger consists of three similar sections, and the internal cavities are connected. The advantage of this structure is that under ideal conditions, each chamber under the same air pressure has the same deformation; thereby, the deformation is repeatable. The research of a single structure is representative, so the deformation simulation of a single cavity is carried out. The constitutive model of the silicone rubber material here adopts the above-mentioned improved Yeoh’s constitutive model. The parameters are detailed in the previous section, and the deformation effect of a single cavity after being loaded by air pressure is programmed to be calculated. Some fixed constraints are used in the upper and lower ends of a single cavity, and air pressure is applied inside. The pressure increases linearly with time. The undamped dynamic equation is computed and simulated from 0 to 0.5 s. Then, the pressure is increased to 12 kPa at the end of 0.5 s.

Figure 13 shows the displacement of the cavity in the Z-axis direction. From the comparison of the results, it can be seen that the expansion effect of a single cavity, calculated based on ANCF and the improved Yeoh model, is good, and the original shape of the cavity is symmetrical. After evenly distributing the air pressure, the deformed structure presents a symmetrical shape, which meets the expected expectation and verifies the correctness of the deformation of a single cavity after being fixedly restrained.

4.3. Analysis of the Overall Deformation of a Flexible Finger

At the top of the finger is a solid encapsulation; during the process of applying air pressure, only the displacement of the top solid section emerges, rather than the deformation. Hence, the air pressure at the top is regarded as an equivalent force applied to the top of the middle of the flexible cavity. The rigid skeleton between the two cavities is simulated by a material with a large elastic modulus. The bottom edge of the finger is fixedly constrained. The calculation method is the same as the above calculation of the single-section cavity. The air pressure is increased linearly with time. The computing time is 0.5 s. At the end of the 0.5 s, the air pressure increases to 12 kPa. The improved Yeoh model is used as the constitutive model. When the air pressure is 8, 10, or 12 kPa, the displacements of the entire finger on the Z-axis are as shown in Figure 14, Figure 15 and Figure 16, respectively.

It can be seen that the overall structure of the flexible finger is effectively bent and deformed under the pneumatic drive. When the air pressure is low, the expansion and deformation, as well as the bending deformation, are not obvious. For instance, the Z-axis expansion displacement is only 7 mm under a pressure of 8 kPa. Additionally, there is nearly no bending. The Z-axis expansion displacement of the finger is increased to 8.2 mm with 10 kPa of pressure. The bending performance can be clearly seen. When the air pressure is up to 12 kPa, the maximum Z-axis expansion displacement is 10 mm. At this moment, the bending angle is larger than at any time. In addition, the expansion and deformation of the flexible finger are symmetrical, which agrees with the predicted tendency. In conclusion, the motion performance of the flexible finger can be characterized by the proposed model effectively, including the expansion, deformation, bending status, and so on, under different air pressures.

5. Conclusions

The ANCF method is known to be well suited for the large deformation analysis of flexible bodies such as flexible fingers. This work presents an improved ANCF low-order plate element based on the hyper-elastic constitutive model, whereby the bending behavior of a silicon rubber flexible finger can be simulated in detail. Compared with the ANCF beam element, conventional lower-order plate element, and higher-order plate element, the novel plate element in this paper characterizes the external contour of the finger better, reflects the bending behavior more realistically, and converges in less computing time.

The effects of volume lock and length constraint were also considered. The former is used to eliminate the incorrect rigidity characteristics of bending deformation for a flexible finger, which is another important factor in realizing large deformation analysis. The skeleton of the flexible finger is simulated by the method of length constraint, which is easy to implement and characterizes the skeleton of the finger.

As a direction of future research, it is desirable to apply the presented method to other ANCF elements, such as solid elements. Furthermore, the method to eliminate the effect of volume lock also needs to be further researched.