An improved dynamic model for a silicone material beam with large deformation

Qiping Xu·Jinyang Liu

Abstract Dynamic modeling for incompressible hyperelastic materials with large deformation is an important issue in biomimetic applications.The previously proposed lower-order fully parameterized absolute nodal coordinate formulation(ANCF)beam element employs cubic interpolation in the longitudinal direction and linear interpolation in the transverse direction, whereas it cannot accurately describe the large bending deformation. On this account, a novel modeling method for studying the dynamic behavior of nonlinear materials is proposed in this paper.In this formulation,a higher-order beam element characterized by quadratic interpolation in the transverse directions is used in this investigation.Based on the Yeoh model and volumetric energy penalty function,the nonlinear elastic force matrices are derived within the ANCF framework.The feasibility and availability of the Yeoh model are verified through static experiment of nonlinear incompressible materials.Furthermore,dynamic simulation of a silicone cantilever beam under the gravity force is implemented to validate the superiority of the higher-order beam element.The simulation results obtained based on the Yeoh model by employing three different ANCF beam elements are compared with the result achieved from a commercial finite element package as the reference result.It is found that the results acquired utilizing a higher-order beam element are in good agreement with the reference results,while the results obtained using a lower-order beam element are different from the reference results.In addition,the stiffening problem caused by volumetric locking can be resolved effectively by applying a higher-order beam element.It is concluded that the proposed higher-order beam element formulation has satisfying accuracy in simulating dynamic motion process of the silicone beam.

Keywords Dynamic modeling ·Nonlinear incompressible hyperelastic materials·Novel modeling method ·Yeoh model·Absolute nodal coordinate formulation

1 Introduction

Flexible material is increasingly attracting attention due to its widespread application in the biomimetic field and industrial engineering.It is crucial to model incompressible materials, which belong to the flexible materials domain,reliably and effectively.For example, bioelastomers,heart muscle and silicone materials can be treated as incompressible materials[1].However,it is necessary to consider both the geometric nonlinearity and the material nonlinearity in modeling nonlinear incompressible hyperelastic material compared with the conventional metal material.Nevertheless,most finite element models adopt the linear constitutive law to study dynamic characteristics of linear elastic materials based on classical theories,which cannot be used in a large deformation event whose internal deformation leads to highly nonlinear dynamic performance similarly to silicone material.An alternative approach to deal with this problem is absolute nodal coordinate formulation(ANCF),proposed initially by Shabana[2],which can be fit for efficient dynamic analysis of flexible material.

The ANCF is a finite element approach in which beam elements can be described with an absolute position and its gradients.In this formulation,gradient coordinates that are partial derivatives of the position vector are used to describe complicated shapes including large rotations and deformations only using a couple of elements.In addition,the ANCF may adopt nonlinear constitutive formulations to model a large amount of nonlinear incompressible material.The ANCF can be used readily with the non-incremental formulation[3].These characteristics makes the ANCF method well suited for the dynamic simulations of the slender, flexible bodies based on nonlinear constitutive models.The use of ANCF leads to a resultant constant mass matrix,and as a consequence,the coriolis and centrifugal forces are identically equal to zero even in the case of spinning structures,which simplifies the description of the equations of motion[4].

Nevertheless,most of investigations devoted to the ANCF assume linear constitutive relation between stresses and strains[5].The nonlinear hyperelastic incompressible material models were introduced into ANCF by Maqueda and Shabana[6]who developed new elastic force formulations of nonlinear materials based on three different nonlinear constitutive laws.These elastic force models can capture some Poisson modes that cannot be captured using other existing beam models in the large deformation analysis.The applications of the incompressible materials to the rubber chains built with ANCF beam elements were presented[7].

ANCF elements that contain all gradient vectors at a node are often referred to as fully parameterized elements.In lower-order ANCF elements[8],several gradients are omitted.However,within fully parameterized elements,different types of locking phenomena may occur due to lower-order interpolation in the transverse direction. In order to overcome this problem, alternative approaches are introduced to formulate the elastic forces[9,10].However,it has been pointed out that fully parameterized beam elements based on general continuum mechanics suffer from Poisson locking[11].Furthermore,Jung et al.[12]suggested and examined the rubber-like material beam using a lower-order ANCF beam element based on three nonlinear constitutive models without highlighting the importance of the volumetric locking phenomena and volumetric locking alleviation techniques.In order to describe the large cross section deformation,a higher-order beam element(HOBE),which has higher-order derivatives as degrees of freedom is introduced.

On the other hand,the fully parameterized ANCF finite elements can also be used to capture the coupling between the cross section deformation and the bending deformation of the finite element[13].However,the cross section deformation modes hinge on the used order of interpolation.The previously developed lower-order beam element(LOBE)with 12 nodal coordinates at each node adopts cubic interpolation in the axial direction and linear interpolation in the transverse direction,which cannot eliminate volumetric locking influence without any specific techniques to deal with volumetric locking phenomenon [14,15]. There are several serious issues including shear locking and volumetric Poisson locking associated with the lower-order interpolation in the transverse direction when the continuum mechanics approach is used[16].

It can be seen that the volumetric Poisson locking phenomenon is of particular interest and can lead to convergence problems regardless of the used number of elements.The previously proposed lower-order fully parameterized ANCF beam element without any specific measures cannot overcome volumetric locking problem.Therefore,it is necessary and important to propose some counter measures to deal with the volumetric locking problem.

This investigation is organized in the following manner.Firstly,an improved dynamic model for nonlinear incompressible hyperelastic materials is developed by employing a higher-order ANCF beam element in order to overcome the volumetric locking problem. Then, the nonlinear elastic force matrices are derived based on the Yeoh model and volumetric energy penalty function. The feasibility and availability of the Yeoh model are verified through static experiment of nonlinear incompressible material. Finally, the effect of adopting three different ANCF beam elements on dynamic response of incompressible hyperelastic material is discussed.Dynamic simulation of the silicone cantilever beam adopting three beam elements is implemented to validate the superiority of the higher-order beam element.It is concluded that the proposed higher-order beam element formulation has satisfying accuracy in simulating the dynamic motion process of the silicone beam.

2 Equations of motion of a three-dimensional beam element

Figure 1 shows a three-dimensional ANCF beam element with two nodes in which shear deformation occurs.In Fig.1,X-Y-Z is the global coordinate system,and x-y-z is local coordinate system,the directions of each coordinate axis in two coordinate systems are shown in Fig.1.

Based on the global shape function and absolute coordinates,the global position vector r of an arbitrary point on the beam can be expressed as

Fig.1 Schematic diagram of three-dimensional beam element

2.1 Higher-order beam element formulation

The existing lower-order ANCF beam element with 24 nodal coordinates is a standard three-dimensional fully parameterized ANCF beam element with two nodes.The absolute coordinate of node k of the LOBE is defined as

where ?rk/?x,?rk/?y,?rk/?z are the position vector gradients of node k in x,y,z directions,respectively.A cubic interpolation in the longitudinal direction and linear interpolations in the transverse directions are adopted to define the LOBE.Nevertheless,the linear interpolation in the transverse direction is insufficient to describe the bending deformation correctly as the volumetric Poisson locking phenomena,which may lead to wrong results.

In order to solve the volumetric Poisson locking problem,the higher-order ANCF beam element with 42 nodal coordinates[17]is adopted for modeling a silicone beam.The absolute coordinate of node k is defined as

where?2rk/?y2,?2rk/?z2,?2rk/?y?z are the second-order derivative vectors in y and z directions,respectively,and the corresponding global shape function matrix is defined as

where I is a 3×3 identity matrix and respective component si(i=1,2,…,14)can be written as

where ξ=x/l,η =y/l,ζ=z/l are the element dimensionless coordinates,and l is the undeformed length of the three-dimensional beam element.

In the present higher-order ANCF beam element formulation,a cubic interpolation is used in the longitudinal direction and quadratic interpolations are used in the transverse directions.Compared with the existing LOBE,the breakthrough in the HOBE is the use of the quadratic interpolation in the transverse directions.The higher-order interpolation can be helpful to avoid the volumetric Poisson locking phenomena and depict bending deformation correctly.

It is important to point out that the shape functions associated with position and gradient nodal coordinates are identical for both LOBE and HOBE.

2.2 Nonlinear elastic force based on Yeoh model

Silicone is a new-type incompressible hyperelastic material,which has a nonlinear stress–strain relationship.In this paper,the strain energy density function is used to describe the isotropic,nonlinear and hyperelastic characteristic of silicone.

Defining J as the matrix of position vector gradients,which is given by

and the right Cauchy–Green deformation tensor C can be written as

The invariants of the right Cauchy–Green deformation tensor can be respectively written as[18]

where tr(C)is the trace of C,det(C)is the determinant of C,and λi(i=1,2,3)are the principle stretch ratios.For an incompressible materials,det(C)=1,or the determinant of the matrix of position vector gradients det(J)=1.

The original model proposed by Yeoh[19]has a cubic form with only I1dependence and is applicable to purely incompressible materials.The Yeoh model for hyperelastic material is a phenomenological model used for analyzing the large deformation of nonlinear incompressible materials such as rubber.The strain energy density function of the Yeoh model is given by Selvadurai[20]

The total strain energy of the nonlinear incompressible Yeoh model is given by

where J=det(J),andκ is penalty coefficient,which should be selected sufficiently large to satisfy incompressibility condition J=1.

Defining rx= ?r/?x,ry= ?r/?y,rz= ?r/?z,from Eq.(6),the determ inant of J can also be written as

In the following section,the nonlinear elastic force formulations of the nonlinear Yeoh model will be proposed to predict the behavior of incompressible hyperelastic materials.Substituting Eqs.(1),(6),and(7)into(8),I1can be written as

where

The virtual work of the elastic force can be deduced by variation of Eq.(12)as

where Qekrepresents the generalized elastic force vector,which can be obtained as

Substituting Eq.(1)into(19),one obtains

In addition,the material constants C10,C20,C30of the Yeoh model can be determined through the uniaxial tensile test.Because of the incompressibility of the nonlinear materials,the stretch ratios have the following relationship

By substituting Eqs.(8)and(21)into(11),the strain energy density function of the Yeoh model also can be expressed in the term of λ as

The relationship between the stretch ratio and the nominal stress of the Yeoh model can be derived by differentiating Eq.(22)with respect to the stretch ratio λ as

2.3 Dynamic equations of motion

On the basis of structural mechanics,the kinetic energy of the beam element can be expressed as

where Meis the constant symmetric element mass matrix,which is defined as

where ρ is the mass density and V is volume of the element in the reference configuration.The coriolis force,centrifugal force,and tangential force as well as other forces are equal to zero in the ANCF.Therefore,the nonzero terms in the dynamic equations of motion derive from the elastic forces and the external forces.

The virtual work of the external force is given by

where Qeerepresents the vector of the element external force.

For a constrained flexible beam with constraint equations Φ=0,based on the principle of virtual work,the dynamic equations of motion of the flexible beam can be expressed as

where M is the global mass matrix,Qkand Qerepresent the vectors of global elastic force and external force,respectively,Φqis the Jacobian matrix,and λ is the vector of Lagrange multiplier.The mass and force matrices are given by

where Beis a Boolean matrix,and the relation between the element coordinate vector qeand the global coordinate vector q is given by qe=Beq.

In order to solve second-order differential Eq.(27)and the constraint equations Φ=0,the New mark numerical integration method combined with the New ton–Raphson method are adopted.Furthermore,generally the Gaussian quadrature is used to integrate over the cuboid-shaped element volume in the finite element analysis[21].

Table1 Material constants of the Yeoh model determined through the uniaxial tensile test

Fig.2 Dimensions of the specimen under the condition of uniaxial tension test(unit:mm)

3 Finite element solutions and analytical solutions under uniaxial tensile tests

First of all,the material constants of the Yeoh model can be extracted from uniaxial tensile tests,and are shown in Table 1.

In this study,the uniaxial tensile test is carried out by using universal testing machine(UTM,Zwick/Roell Z050)to determine the material constants of the Yeoh model by minimizing the error between the stretch ratio and the nominal stress based on the principle of the least squares method.The specimen is tested for a couple of times,then the obtained strain–stress curve is converted into the nominal stress–stretch ratio curve.

Figure 2 illustrates the dimension of the specimen under the condition of uniaxial tension test.As shown in Fig.2,W,T,L,and R are width,thickness,gauge length,and radius of shoulder of the specimen,respectively.A tensile speed of 50 mm/min is applied to the specimen,and the test ends when the stretch ratio reaches 4.2.Finally,the experiment data which is denoted by a small blue circle is shown in Fig.3b.

The elastic forces determined by Eq.(17)are calculated by using the material constants.Figure 3 shows the stress–stretch ratio curves gained by comparing the results obtained using the Yeoh model with the test data.A horizontal static load is applied to the right end of the silicone beam in the step-by-step loading manner.The equilibrium equations combined with the constraint equations are given by

Fig.3 The stress–stretch ratio relationship curves.a Yeoh model,b Yeoh model and experiment data

The nonlinear Eq.(29)is solved by using New ton–Raphson iteration technique.Finite element solutions obtained through solving nonlinear Eq.(29)are compared with analytical solutions acquired solving Eq.(23).Stress–stretch ratio relationship curves of the Yeoh model are shown in Fig.3a,and the stress–stretch ratio relationship curves of finite element solutions for the Yeoh model and experiment data are shown in Fig.3b.

As can be seen from Fig.3a, finite element solutions are nearly consistent with analytical solutions for the Yeoh model in the case of the number of elements is enough(such as 60 elements),and the accuracy of the finite element model is verified.Moreover,Fig.3b is the stress–stretch ratio curves of finite element solutions for the Yeoh model compared with the experiment data,and it is revealed that the Yeoh model can well capture the experimental data.

Fig.4 The silicone cantilever beam model in the absolute nodal coordinate formulation(unit:mm)

Fig.5 The silicone cantilever beam in the static equilibrium

4 Static experiment and dynamic simulation

4.1 Static experiment

Firstly,static equilibrium experiment of the silicone cantilever beam is implemented.Figure 4 shows the silicone cantilever beam model.The initial length of the beam(L)is 160 mm and the width(b)and height(h)of the cross section are7and 5mm,respectively.The mass density of the beam is 1100kg/m3.The material constants shown in Table1are used to calculate the elastic force of Eq.(17).The incompressibility constant used in Eq.(12)is assumed to be κ=1000 MPa considering the incompressibility.

Figure 5 shows the silicone cantilever beam which is in the static equilibrium state only under the gravity force. The silicone beam is clamped at the reference plane. Mark points are located at the points 50, 100, 150, and 160mm from the reference as shown in Fig. 6. Then, static simulation experiment of the silicone cantilever beam through MATLAB mathematical software and ANSYS engineering software is performed.

Fig.6 Marks points on the silicone beam

Fig.7 The static equilibrium state of the silicone cantilever beam

Figure 7 shows the reference result obtained using the ANSYSBEAM 188element and the solution acquired using the HOBE and experimental data for the silicone cantilever beam,which is in static equilibrium only under the force of gravity.The red solid line is the ANSYSBEAM 188element result,the blue dashed line is the HOBE solution,and the green hexagram is the experimental data.As can be seen from the Fig.7,the ANSYS BEAM 188 element result matches the HOBE solution well.However,the experimental data lag behind the ANSYS BEAM 188 element result and the HOBE solution,possibly because of the measurement error.It is disclosed that there is a good agreement between the HOBE solution and the ANSYS BEAM 188 element result as reference result overall,and both the HOBE solution and the ANSYSBEAM 188element result can capture the experimental data accurately.

In addition,the dynamic simulation of the silicone cantilever beam only under the action of gravity force is carried out using the ANSYS BEAM 188 element and adopting LOBE and HOBE,respectively.

4.2 Dynamic simulation

The equations of motion and the constraint equations of a constrained flexible beam are given by

As previously mentioned,herein,the LOBE with 24 nodal coordinates refers to the standard three-dimensional fully parameterized ANCF beam element without any techniques to solve the volumetric locking issue.The LOBE and the HOBE are employed to investigate the dynamic motion process of the silicone cantilever beam.The reference solution can be obtained using the ANSYS BEAM 188 element.This BEAM 188 element have six or seven degrees of freedom per node,the first six degrees of freedom account for shear deformations and the cross-section deformations,while the seventh degree of freedom allows for warping.The BEAM 188element is based on Timoshenko Beam theory.At the clamped end,all nodal coordinates are fixed.Nevertheless,the additional degrees of freedom related to higher-order terms at the clamped end can be released without observably affecting the results[22,23].

Figure 8 shows the configurations of the silicone cantilever beam at two different moments acquired from LOBE,HOBE,and the ANSYS BEAM 188 element based on the Yeoh model.The red hexagram is the ANSYSBEAM 188element model result,the green solid line is the HOBE result,and the blue dashed line is the LOBE result.Figure 8 is the comparative configurations gained using the LOBE,the HOBE,and the ANSYS BEAM 188 element.As can be seen,the result acquired utilizing the HOBE accurately coincides with the result gained with the ANSYS BEAM 188 element as reference result overall.Both results are in a good agreement,whereas the result acquired using the LOBE is inconsistent with the reference result.It has been proved that adopting the HOBE is valid for the dynamic characteristics analysis and actual dynamic motion process of the silicone cantilever beam,while employing the LOBE is invalid.

Figure9 shows the variations of the horizontal and vertical displacements of the mark point located at x=150 mm over time obtained employing the LOBE based on the Yeoh model.As the increase of the deformation,the horizontal and vertical displacements of the mark point located at x=150 mm fluctuates within a certain range(i.e.1.4 mm),and it is shown that high-frequency oscillation occurs within a small amplitude.The reason for this result may be attributed to the volumetric locking which occur in nonlinear incompressible hyperelastic materials.

Fig.8 The configurations of the silicone cantilever beam.a t=0.1 s,b t=0.5 s

Figure 10 shows the variations of the horizontal and vertical displacements of the mark point located at x=150 mm over time achieved from the LOBE,HOBE,and the ANSYS BEAM 188element.There d hexagram is the ANSYSBEAM 188 element result,the green solid line is the HOBE result,and the blue dashed line is the LOBE result.As can be seen from Fig.10,the horizontal and vertical displacements acquired employing the HOBE can substantially follow the result gained by the ANSYSBEAM 188element as reference result overall,both results are in a good agreement,whereas the results gained using the LOBE thoroughly fails to follow the reference result.It also has been validated that adopting the HOBE is valid and workable for actual dynamic motion process of the silicone cantilever beam,while employing the LOBE is invalid.The use of the HOBE in simulating dynamic motion process of the silicone beam suppresses volumetric locking and achieves satisfactory results,whereas the stiffening problem caused by volumetric locking cannot be solved by the LOBE.It is further suggested that the Yeoh model has satisfying accuracy in simulating actual dynamic motion processes of the silicone beam.

Fig.9 Horizontal and vertical displacements of the mark point located at x=150 mm.a Horizontal displacement,b vertical displacement

4.3 Comparison with modified lower-order beam element

The modified lower-order beam element(MLOBE)[24]is the commonly used beam element with selective reduced integration approach. For this method, the volumetric energy penalty function item, which allows for the volumetric behavior is calculated only along the beam centerline without considering integration in y and z directions. Figure11shows the difference of the HOBE and MLOBE at t=0.1s.The red hexagram is the ANSYSBEAM 188 element result,the green solid line is the HOBE result,and the violet dashed line is the MLOBE result.It can be seen that the results obtained by HOBE are more consistent with the results acquired by the BEAM 188 element than the MLOBE.

Fig.10 Comparison of horizontal and vertical displacements of the mark point located at x=150 mm.a Horizontal displacement,b vertical displacement

Fig.11 Comparison of two beam elements

5 Conclusions

In this study,an improved dynamic model for nonlinear incompressible hyperelastic materials is proposed by using HOBE characterized by cubic interpolation in the longitudinal direction and quadratic interpolation in the transverse directions.This model can describe the large bending deformation effectively based on the Yeoh model and volumetric energy penalty function.The elastic force formulations of the Yeoh model with nonlinear constitutive law are derived considering the large deformation of flexible material utilizing ANCF beam element.Dynamic simulation of the silicone cantilever beam only under the action of gravity adopting the LOBE model and HOBE model is implemented.It is noteworthy that the proposed HOBE formulation has satisfying accuracy in simulating the dynamic motion process of the silicone beam.

The effect of employing three different three-dimensional ANCF beam elements,namely LOBE,MLOBE and HOBE,on the dynamic response of silicone material in the incompressible domain is investigated.It is shown that exploiting HOBE is crucial for executing reliable dynamic analysis of incompressible material based on the nonlinear constitutive model.The simulation results obtained based on the Yeoh model applying two existing ANCF beam element models are compared with the results achieved from a commercial finite element package as reference results.It is found that the results acquired using the HOBE model have good agreement with reference results,while the results obtained utilizing LOBE are entirely different from the reference results.The stiffening problem caused by volumetric locking can be solved effectively by adopting HOBE.Furthermore,it is demonstrated that HOBE is more accurate than MLOBE.

It was also validated that employing the HOBE model is valid and workable for the actual motion process of silicone cantilever beam.In addition,the static equilibrium experiment of the silicone cantilever beam under gravity is carried out.It has been proved that the feasibility and availability of the Yeoh model are verified in the large deformation analysis for nonlinear materials.It is further revealed that the Yeoh model has satisfying accuracy in simulating the dynamic motion process of the silicone beam.

AcknowledgementsThis research was supported by the National Natural Science Foundation of China(11772186 and 11272203).