Large deformations of hyperelastic curved beams based on the absolute nodal coordinate formulation

Compared with traditional linear elastic materials, the soft structure composed of incompressible hyperelastic materials has not only geometrical nonlinearity but also material nonlinearity during deformation. In this paper, the absolute nodal coordinate formulation (ANCF) is used to study the large deformations and large overall motions of incompressible hyperelastic curved beams. A novel large deformation dynamic modeling method for curved beams made of hyperelastic materials is proposed, in which a simplified Neo-Hookean model is combined with the one-dimensional ANCF beam element. The elastic force vector is calculated according to the exact expression of curvature. The dynamic equations are derived by using the virtual work principle. The dynamic responses of a cantilever silica gel beam under gravity are calculated based on the present method and compared with those of the improved low-order beam element (ILOBE), high-order beam element (HOBE), and commercial finite element analysis software (ANSYS). Simulation results show that the proposed method can accurately describe the large deformation and large overall motion of the beam, and has better computational efficiency. Research in this paper provides an efficient dynamic model for the dynamics analysis of soft robot arms.


Introduction
Soft robots are usually made of hyperelastic soft materials or incompressible materials (e.g., rubber, silicone) that can be continuously deformed (similar to muscle drive) [1], which can achieve large deformations and large overall motions [2]. In order to more accurately describe the dynamical behavior of soft robots, it is important to consider not only the geometric nonlinearity due to the large range of motion but also the material nonlinearity due to the application of hyperelastic materials. Therefore, it is a hotspot in recent years to establish an accurate dynamics model for soft body structures using multibody dynamics theory and to study efficient and stable methods for calculating the dynamical responses.
In the past decades, abundant research results have been achieved in this field, and many practical engineering problems have been solved. A large number of practical engineering components (such as flexible robot arms, flexible satellite antennas, large span flexible connecting rods, etc.) can be simplified into the mechanical model of flexible beam [3]. Due to the influence of new materials and technologies, the large deformation of the mechanism becomes more obvious. Shabana [4] proposed the absolute nodal coordinate formulation (ANCF) in 1996. The ANCF describes the motion of a flexible body in a unified inertial coordinate system and avoids the complex cornering problem of the flexible body using the slope vector. There will be no Coriolis force and centrifugal force, and the elastic force is the only nonlinear term. The simplicity of programming in solving large deformation problems has greater advantages and has achieved more widespread applications [5,6].
Modeling for curved beams can be broadly divided into two types, one using straight beams to approximate curved beams and the other developing new curved beam elements. Shabana [6] first applied ANCF to flexible bodies with discontinuous slope vector coordinates and demonstrated that the method can effectively analyze large deformations of flexible curved beam structures through the example of a folded beam. Since the elements established by Shabana may produce coupled deformation and shear locking, Sugiyama et al. [7,8] adopted the Hellinger-Reissner variational principle to correct the shear stress and established a fully parametric curved beam element and gradient defect curved beam element for analyzing large deformations in flexible multibody systems. Zhang et al. [9] proposed Bernoulli-Euler beam elements with high-precision curvature constraints based on ANCF, which can analyze nonlinear problems of large deformation and large rotation of planar thin beam structures. Besides, Guo et al. [10] and Wu et al. [11] used ANCF to derive a onedimensional two-node curved beam element based on the exact expression of curvature and studied the large deformation problem. Shabana et al. [12] used strain splitting method to study the Poisson locking problem of curved beam and plate structures, and compared the results with commercial finite element software to prove the validity. Based on the work of [11], Hewlett et al. [13] developed a more efficient first-order integration method. It should be mentioned that, all the above studies are based on linearly elastic materials.
Most current curved beam models are limited to the linear constitutive equation with stress-strain, however, hyperelastic materials usually have nonlinear constitutive equation and are usually completely incompressible [14,15], and all commercial finite element method (FEM) packages for nearly incompressible materials do not fully satisfy this constraint [16].
Farzam et al. [16,17] developed a fully Lagrangian finite element calculation method based on the threedimensional constitutive equations of isotropic hyperelastic materials. However, their study used traditional generalized coordinates, which made it difficult to capture some important modes in large deformations. Maqueda et al. [18] and Jung et al. [19] first used the ANCF to derive the elastic force expressions for three nonlinear hyperelastic material constitutive models (Neo-Hookean, Mooney-Rivlin, and Yeoh), and applied them to three-dimensional low-order beam element and discussed the validity and simulation accuracy of the three models through simulation and experiment. Luo et al. [20] established hyperelastic thin-shell elements of compressible Neo-Hookean and incompressible Mooney-Rivlin materials and compared the shell dynamics response of different materials. Orzechowski and Frączek [21,22] based on their work investigated a volume-locking suppression method applied to a nonlinear hyperelastic nearly incompressible material model with fully parametric beam elements. Based on the work of [21,22], Xu et al. [23,24] improved the low-order beam element (ILOBE) using a reduced integration technique and proposed an ANCF high-order beam element (HOBE) with quadratic interpolation in the transverse direction. Two kinds of beam elements are combined with four kinds of nonlinear material models, and the accuracy of the different nonlinear material models is compared by numerical simulation arithmetic and physical experiments. In addition, Greco et al. [25] comprehensively compared the classical strain measurements used in the hyperelastic material model and proposed a new Generalized Hyperbolic Sine (GHS) strain measurement.
However, in the previous work, researchers all used three-dimensional fully parametric beam elements to simulate the bending large deformation of straight rubber beams, and the derivation of the elastic force vector was also based on the three-dimensional nonlinear constitutive equation of continuum mechanics. The elements had more degrees of freedom and the nonlinear material itself requires more iterations to converge, which leads to computational inefficiency.
In this paper, a new dynamical modeling method for large deformations and large overall motions of hyperelastic material curved beams is proposed. The ANCF one-dimensional two-node curved beam elements are combined with the Neo-Hookean model simplified to the uniaxial tensile case, and the elastic force vector with precise curvature expression and nonlinear constitutive relation is derived. In Sect. 2, the dynamic equation based on ANCF is established. The internal force of the system is nonlinear elastic force. In Sect. 3, four numerical examples of statics and dynamics are given to verify the correctness and effectiveness of the proposed method. Compared with other methods, this method can not only describe the large deformation of hyperelastic material beams with initial curvature, but also has higher efficiency, and has certain applicability to different nonlinear constitutive relations.

Curved beam element based on ANCF
The large deformation a flexible Euler Bernoulli curved beam is shown in Fig. 1. The initial central angle of the curved beam is a 0 and the initial radius of curvature is R s , which is much larger than the height of the cross section. The cross section of the beam is perpendicular and symmetric to the central axis, and the rotational inertia and shear deformation of the beam cross section is not considered. The flexible pendulum beam is discretized by finite elements, and the beam with length L is equally divided into n curved beam elements. Figure 2 shows the schematic diagram of the curved beam element before and after the deformation, where OXY is the global inertial coordinate system and O e X e Y e is the element coordinate system.
The global absolute position vector r of any point P on the axis of the curved beam element can be expressed as: where r 1 ; r 2 is the component of the absolute position vector of the point P in the X and Y directions, respectively; and x is the arc coordinates of the point P on the center axis of the curved beam element in the element coordinate system. The matrix S is the shape function of the curved beam element [11,12] as: in which where n ¼ x=l, and l is the length of the center axis of the beam element without deformation in Fig. 2.
The absolute nodal coordinate vector e of the 8-DOF curved beam element in the inertial coordinate system can be expressed as:

Nonlinear elastic forces based on the Neo-Hookean model
Hyperelastic materials are distinguished from linear elastic materials by no longer following the small deformation assumption and having nonlinear constitutive relations. The deformation energy of a material can be written as a scalar function related only to its current strain state. The strain energy density function U of a material can be expressed as a function of three invariants I 1 ; I 2 ; I 3 that depend only on the right Cauchy-Green deformation tensor C as [26]: where the three invariants of the right Cauchy-Green deformation tensor C can be expressed as: where k i ¼ 1 þ e i is the tensile ratio in each of the three principal axes and e i is the strain in the tensile directions. The partial derivative of strain energy density function U with respect to the principal tensile ratio k i can obtain the relationship between stress and strain of hyperelastic material as [27]: The strain energy density function of the oneparameter Neo-Hookean model for hyperelastic materials can be expressed as [18,28]: where C NH is the material coefficient related to the deformation response, which is usually measured experimentally, k 0 is the volume modulus of the material to ensure the incompressibility of the material, and J is the determinant of the deformation According to the uniaxial tensile test characteristics, the relationship between the principal tensile ratio is: By substituting Eqs. (9) and (8) into Eq. (7) and assuming that the material is completely incompressible, the nonlinear stress-strain relationship in the tensile direction of the Neo-Hookean model is obtained as According to continuum mechanics, the strain e at any point on the non-central axis of the curved beam element can be expressed as the expression of the exact curvature of the curved beam element as: where e 0 represents the longitudinal strain on the axis of curvature, which can be replaced by the Green Lagrange strain: j is the curvature change before and after deformation at the centerline, which can be expressed by the partial derivative of the element's absolute position vector r about the curved beam arc coordinate x: in which Substituting Eq. (1) into Eq. (12) and Eq. (13) can be obtained in which Based on the Euler-Bernoulli beam hypothesis, the virtual work of the elastic force of the curved beam element based on the Neo-Hookean model of the incompressible hyperelastic material is: in which V is the volume of the beam element. Substitute Eq. (11) into Eq. (18), and assume that the hyperelastic beam is symmetric about the central axis, then R y dA ¼ 0, and the moment of inertia of the cross section is I z ¼ R A y 2 dA. The generalized elastic force of the element can be expressed as: It is assumed that B is the Boolean positioning matrix in the beam element absolute node coordinate vector e corresponding to the total absolute node coordinate vector q, e ¼ Bq can be obtained. Then the virtual work done by the total elastic force of the beam system is: where Q k is the generalized elastic force vector of the beam system as:

Generalized mass matrix and generalized external force vector
The kinetic energy of the beam element can be expressed as: where M e is the beam element mass matrix, which can be expressed as: in which q 0 is the density and A is the cross-sectional area of the beam element without deformation. Thus, the total kinetic energy of the beam can be expressed as: where the overall generalized mass matrix M of the system is: Beam elements are subjected to uniformly distributed gravity, and the virtual work done by gravity can be expressed as: The curved beam moves only in the vertical plane and has G ¼ ½0; Àqg T in the inertial coordinate system, where g is the gravitational acceleration; Q ek is the gravity vector of the element, which can be expressed as: So the generalized virtual work done by the external force of the curved beam is: The generalized external force vector Q g of beam system is

Dynamic equations
By adopting the virtual work principle and combining Eqs. (20), (24) and (28), yields: The boundary conditions of the curved beam are considered and transformed into the constraint equation U q; t ð Þ ¼ 0. By the Lagrange multiplier method, the set of differential-algebraic dynamics equations of the flexible curved beam multibody system can be obtained: When the statics problem is solved, the statics equilibrium equation becomes:

Numerical simulations
In order to verify the feasibility and correctness of the proposed method in analyzing the large deformation problems of incompressible hyperelastic curved beams, statics and dynamics simulations of flexible curved beams subjected only to gravity are carried out through several examples in this section. In numerical calculation, the Gaussian integral method was used to integrate the volume of the curved beam element, the Generalized-a method [29,30] was used to solve the dynamics Eq. (31), and the Newton-Raphson method was used to solve the static Eq. (32).

Statics simulation of a cantilever curved beam
Consider the cantilever silica gel curved beam model as shown in Fig. 3. The left end of the curved beam is fixed and the right end of the curved beam is free. The curved beam falls freely from the initial position shown in the vertical plane. The physical parameters of the beam are shown in Table 1.
The number of elements affects the convergence of the ANCF calculation results. Different numbers of elements are selected to calculate the same problem. Figure 4 shows the absolute position coordinates of the tip of the curved beam (the right node of the last element) in the horizontal and vertical directions during the static equilibrium, respectively. When the number of elements is less than 10, the numerical error cannot be ignored, and when the number of elements reaches 12, the simulation results of the free end of the curved beam are sufficiently converged, which also shows that the method proposed in this paper does not have locking problem.
To further verify the correctness of the method in this paper, the above model is discretized into 14 elements to simulate the configuration of the silica gel curved beams in static equilibrium, and the simulation results are compared with those of the commercial finite element software ANSYS APDL. Because the BEAM element in ANSYS is a one-dimensional element, while the hyperelastic material model is based on three-dimensional analysis, using BEAM element ANSYS software will show that the analysis   Fig. 4 Variations of a the horizontal and b the vertical displacements of the tip of the curved beam with different numbers of elements cannot be performed. Therefore, SOLID element is selected in ANSYS. Since the curvature of the hyperelastic material BEAM is significant, SOLID186 element with a higher approximate polynomial order is selected which is suitable for large deformation and large overall motion analysis of the Neo-Hookean constitutive model of hyperelastic materials. The Neo-Hookean material model was selected in ANSYS analysis, and the geometric parameters and material parameters are shown in Table 1. The finite element model of the curved beam is shown in Fig. 0.3 (b). It is divided into 320 finite elements with 80 equal parts along the length and 2 equal parts along the width and height of the curved beam, and large deformation is allowed during simulation. Figure 5 shows the comparison of the simulation calculation results between the proposed method and the ANSYS, where the red solid line is the initial configuration of the curved beam, the black solid line is the numerical calculation result of the proposed method, and the blue-dashed line is the result using ANSYS. Table 2 shows the comparison of the coordinates of the curved beam between ANSYS and current method at different positions (x ¼ L=4, It can be found that the results obtained by the current method are consistent with those calculated by ANSYS. The errors between the beam position coordinates calculated by current method and the position coordinates simulated by ANSYS is small. Therefore, the method proposed in this study can accurately describe the large deformation problem in the static analysis of hyperelastic materials.

Dynamic simulation of a cantilever curved beam
In order to verify the simulation accuracy of the proposed method in analyzing the dynamics of hyperelastic material curved beam, one dynamics simulation is performed for the model shown in Fig. 3.
The initial central angle of the hyperelastic curved beam a 0 is p=3 rad, the radius of curvature R s is 1 m, the section size is 0:05 m Â 0.05 m, and other parameters are consistent with Table 1. In ANSYS analysis, the parameters are consistent, allowing large deformation during simulation. The finite element model of the curved beam is similar to Fig. 3 (b). The curved beam is divided into 1080 finite element elements by using SOLID186 element, with 120 equal parts along the length direction of the curved beam and 3 equal parts along the width and height direction. The simulation is calculated for 1.0 s. Numerical damping is not considered in the generalized-a calculation. Figure 6 shows the evolution of the absolute position of the endpoint over time, where the black and blue solid lines indicate the absolute positions in the horizontal and vertical directions of the proposed method respectively, and the red and green-dashed lines indicate the calculation results of ANSYS, respectively. Figure 7 shows the deformed configurations at 0.2 s interval using the method of this paper compared with the deformation at the corresponding moment in ANSYS. Table 3 shows the comparison of the tip coordinates of the curved beam based on the current method and ANSYS simulation at different times.
It can be found that there is a good consistency between the results, and the simulation results of the proposed method and ANSYS simulation results converge to almost the same solution. With the increase in time, the deformation of the silica gel beam increases, and there is a small error between the two calculation results, but the current method still has a good simulation accuracy. The calculation results of the current method in this paper are consistent with the commercial software ANSYS. It is worth noting that the number of elements divided by the method proposed in this paper is 16, which is much less than Fig. 5 The configuration of the curved beam that of the ANSYS solid model that requires 1080 elements. In general, this method can accurately calculate the dynamic response of beams with large deformation and large overall motions.

Dynamic simulation of the straight cantilever beam
The method in this paper is based on the exact curvature expression. When the initial curvature of the curved beam is j 0 ¼ 0 and the initial node coordinate vector in the form of the straight beam is adopted, the curved beam model degenerates into straight beam model. Dynamic simulation was carried out on the cantilever silica gel straight beam fixed at the left end. The physical parameters of the straight beam are consistent with Ref. [21], where the coefficient of hyperelastic material C NH is 0:8 MPa, the density q is 7200 kg/m 3 , the length of the straight beam L is 1 m, and the section size is 0:02 m Â 0.02 m. In ANSYS simulation, since the structure has no initial curvature, SOLID185 finite element is used to discretize the straight beam into 320 elements with 80 equal parts along the length and 2 equal parts along the width and height of the beam. The material and geometric parameters are consistent with the current method, allowing large deformation. The simulation results of the proposed method are compared with the results calculated by the ILOBE in Ref. [23], the HOBE in Ref. [24], and ANSYS. It should be noted that the elastic force vector of hyperelastic materials modeled by ANCF three-dimensional ILOBE and HOBE are all based on the three-dimensional constitutive model of the right Cauchy-Green deformation tensor. Penalty functions are required to ensure incompressibility, and the penalty coefficient k ¼ 1000 MPa is selected [24]. Figure 8 shows the comparison of the calculation results of the dynamics simulation using four different methods over time, showing the absolute positions in horizontal and vertical directions, respectively. Among them, the black solid line is obtained by using  . 6 The absolute displacements of the tip of the curved beam Fig. 7 The deformed configurations of the curved beam from 0 s to 1.0 s the method proposed in this paper, the red dash-dotted line is obtained by using ILOBE, the green double dots line is obtained by using HOBE, and the blue dashed line is obtained by using ANSYS simulation results. It can be found that there is a good consistency between the results, and all the results can converge to almost the same solution. When t 1:0 s, there is little difference between the results of the four methods, which indicates that the four methods can describe the large deformation of the silica gel straight beam well. However, as time continues to increase, the difference between HOBE and other methods increases. This is because HOBE elements have more degrees of freedom, and more elements are needed to ensure the computational accuracy for large deformation. Table 4 shows the elements required and the CPU runtime consumed by the four different methods using the same computational platform (the operating memory is 8 GB, the processor is Inter Core I5-8300H, and the CPU frequency is 2.3 GHz) to calculate the large deformation of the straight beam for 1.5 s. It can be found that when the model both established based on ANCF is used for dynamic numerical simulation, the proposed method takes the least number of elements and consumes the shortest CPU time to calculate the same problem, which is nearly 1/4 of ILOBE and 1/40 of HOBE. Therefore, the proposed method can not only accurately describe the large deformation and large overall motions of the beam with fewer elements, but also has better computational efficiency compared with other methods.

Dynamic simulation of a curved beam pendulum
The Neo-Hookean model has some limitations under large strains, and more hyperelastic constitutive models are worth considering. In order to verify the correctness of the proposed method when dealing with curved beams with different material constitutive models and different geometric parameters, a silica gel curved beam pendulum with the hinged left end and the free right end was considered as shown in Fig. 1. The curved beam pendulum falls under gravity and two incompressible material models are used for the calculation: Neo-Hookean and Mooney-Rivlin [18,28]. The stress-strain relationship of Mooney-Rivlin is shown in Eq. (33), where the material coefficient is C 1 ¼ 0:8 MPa and C 2 ¼ 0:2 MPa respectively, the density of the hyperelastic curved beam q is 7200 kg/m 3 , the initial central angle a 0 is p=3 rad, the length of the curved beam L is p=3m, and the section size is 0:05 m Â 0.05 m.
Note that when C 2 ¼ 0, this model is degraded into the previously introduced Neo-Hookean material model. In addition, when only small strains are considered, the linear elastic material model can degenerate, and Young's modulus is E ¼ 6ðC 1 þ C 2 Þ [21]. Figure 9 shows the comparison of the absolute position of the tip in horizontal and vertical directions over time with different material models. The black line, blue line, and red line are simulation results based on the Neo-Hookean nonlinear constitutive model, Mooney-Rivlin nonlinear constitutive model, and equivalent linear elastic model under small strain, respectively. It can be seen that when t\0:5 s, the calculation results of the three are not much different. However, with the increase in time, the simulation results based on the Neo-Hookean model and Mooney-Rivlin model are still in good agreement, whose dynamic characteristics are similar, but they are quite different from the linear elastic constitutive model. This is because the small strain hypothesis of linear elasticity is no longer satisfied with the increase in deformation. Therefore, incompressible nonlinear material constitutive relations need to be considered in the analysis of such large deformation problems, which also reflects the applicability of the proposed method to different constitutive models of hyperelastic materials.
In addition, the geometric parameters of the material also have a great influence on the deformation of curvature. For intuitive analysis, it is necessary to transform the absolute node displacement into transverse deformation and longitudinal deformation of the tip. Suppose the rotation angle of the coordinate system O e X e Y e relative to the coordinate system OXY in the curved beam model is u, then: where u 1 ; u 2 is the longitudinal and transverse deformation of the tip of the curved beam model, respectively, and H is the direction cosine matrix: The simulation model is shown in Fig. 1. In order to control the variables, the length of the curved beam L remains unchanged, and the initial curved angles a 0 are set as p=6 rad, p=4 rad, p=3 rad and p=2 rad, respectively. Figure 10 shows the comparison of longitudinal and transverse deformation of the end of a curved beam pendulum at different initial central angles over time. It can be found that within 1.5 s, the maximum longitudinal deformation and transverse deformation of the curved beam pendulum system will increase with the increase in the initial curved angle at any time, and the response time to reach the maximum longitudinal deformation and transverse deformation will shift backward. Therefore, more attention should be paid to the problem of the large deformation of hyperelastic curved beams with large curvature.

Conclusions
In this paper, a new dynamic modeling method for hyperelastic materials with large deformation curved beams is proposed based on the ANCF one-dimensional beam element. The elastic force vector established by this method is based on a simplified nonlinear incompressible Neo-Hookean model and the accurate curvature expression. Several statics and dynamics simulation examples of a cantilever silica gel curved beam under gravity are calculated by using this method, and the correctness and effectiveness of the model are validated. The calculation results based on the proposed method are consistent with those of the improved low-order beam element (ILOBE), highorder beam element (HOBE), and commercial finite element analysis software (ANSYS), and the current method requires fewer elements and is more efficient. In addition, when the deformation is too large, linear elastic constitutive relation is no longer applicable, and nonlinear material constitutive relation needs to be adopted. The method in this paper has certain applicability to different nonlinear constitutive models of hyperelastic materials.
The dynamic model based on ANCF proposed in this paper can accurately and efficiently simulate the dynamics of large deformation and large overall motion of hyperelastic curved beams, and can improve the computational efficiency of hyperelastic curved beams and extends the previous work only on superelastic straight beams to curved beams. However, the above methods and experimental studies still need to be further improved to develop more efficient solution methods to deal with more refined models. It is hoped that this dynamic model can be extended to the optimal design and dynamic control of soft robots and flexible manipulators.