Large Sliding Contact Analysis Of A Buckle

Chapter 27: Large Sliding Analysis of a Buckle

27

Large Sliding Contact Analysis of a Buckle 

Summary



Introduction



Modeling Details



Solution Procedure



Results



Modeling Tips



Input File(s)



Video

463 464 464

469

473

472 473

467

CHAPTER 27 463 Large Sliding Contact Analysis of a Buckle

Summary Title

Chapter 27: Large Sliding Contact Analysis of a Buckle

Features

Deformable-deformable contact, bilinear, Coulomb friction model, Hookean, isotropic elastic material, adaptive time stepping, solid elements with assumed strain formulation

Geometry

168 mm

Y

X Z

247

y metr Sym Half

mm

Material properties

E = 10GPa ,  = 0.4

Analysis characteristics

Quasi-static analysis using: adaptive time stepping and geometric nonlinearity due to large displacement

Boundary conditions

Sliding, frictional contact with: ends fixed for second contact body and contact between the two deformable bodies with  = 0.1

Applied loads

Prescribed displacements for the end nodes of the first contact body with two load cases: insertion (clipping) and removal of the buckle

Element type

8-node solid element with assumed strain formulation

FE results

1. History plot of y-displacements for specific nodes 2. Normal and frictional contact force comparison of Nastran and Marc 3. Load displacement curves comparison between the frictional and frictionless cases Fx

Fx

1000

Fx (N)

500 0.5

0

1.0

1.5 Time (s)

-500

Frictionless Frictional

-1000 -1500

Insert

Remove

-2000 Fx

Fx

464 MD Demonstration Problems CHAPTER 27

Introduction This problem demonstrates the ability of MD Nastran SOL 400 to do a frictional contact problem. An ostensibly simple geometry poses a substantial challenge for the contact algorithm due to the large sliding involved between the two deformable bodies. Sudden changes in the motion path pose a challenge to the ability of the contact algorithm to correctly place the node on the contact surface while respecting the various geometric details in the problem. Due to large bending stresses in the deformed configuration, assumed strain formulation is used with the 8-node hexahedral elements. The material is elastic and the original geometry without residual stresses is recovered upon the complete removal of the loading. From elementary strength of materials analysis, the tip deflection for beam bending can be written as:  = PL 3   3EI 

where P is the applied load, L is the length of the beam, I is the moment of inertia and E is the Young’s modulus. The normal stress along the beam cross section varies in the thickness direction as:  xx = M t  I

where M is the moment and  t is the thickness coordinate. It must be noted that the above solution only holds for small displacements and uniform cross section.

Modeling Details A numerical solution has been obtained with MD Nastran’s SOL 400 for a 3-D representation of a belt buckle with a deformable-to-deformable contact between the two pieces of the buckle. The details of finite element model, contact simulation, material, load, boundary conditions, and solution procedure are discussed below. The case control section of the input contains the following options for nonlinear analysis: SUBCASE 1 STEP 1 TITLE=Insertion (Clipping) ANALYSIS = NLSTATIC NLPARM = 1 BCONTACT = 1 SPC = 2 LOAD = 1 DISPLACEMENT(PLOT,SORT1,REAL)=ALL SPCFORCES(PLOT,SORT1,REAL)=ALL STRESS(PLOT,SORT1,REAL,VONMISES,BILIN)=ALL NLSTRESS(PLOT,SORT1)=ALL STEP 2 TITLE=Removal ANALYSIS = NLSTATIC NLPARM = 2 BCONTACT = 2 SPC = 6 LOAD = 2 DISPLACEMENT(PLOT,SORT1,REAL)=ALL SPCFORCES(PLOT,SORT1,REAL)=ALL STRESS(PLOT,SORT1,REAL,VONMISES,BILIN)=ALL NLSTRESS(PLOT,SORT1)=ALL

CHAPTER 27 465 Large Sliding Contact Analysis of a Buckle

The analysis contains a single subcase with two steps. The two steps comprise of individual load sequences consisting of insertion (clipping) and removal of the belt buckle. Each step has a definition of convergence control option via NLPARM, contact table and parameters via BCONTACT, applied displacements (or single point constraints) via SPC and the displacements and stress results for the .f06 (output) file. A zoomed-in view of the cross section of the model shown in Figure 27-1 consists of an outer piece modeled as body 2, the buckle, while the inner piece is modeled as body 1, the insert.

Figure 27-1

Geometry and a Zoomed-in View of a Belt Buckle

Large displacement effects are included in the nonlinear analysis using the option: PARAM

LGDISP

1

While the assumed strain formulation is flagged using the option: NLMOPTS,ASSM,assumed

The NLMOPTS field triggers the assumed strain formulation which provides a better bending behavior of the continuum elements. This alleviates the difficulty associated with spuriously large shear stresses induced due to bending moment. The LGDISP field indicated the use of large displacement, large rotation kinematics of the element. This is adequate when the analysis consists of Hookean elastic material; however, incase of large deformation plasticity or other inelastic models, the LRGSTRN parameter should be used in the NLMOPTS option (for more details on its usage, please refer to : Chapter 3: 3-D Sheet Metal Forming of this manual).

Element Modeling Besides the standard options to define the element connectivity and grid coordinate location, the bulk data section contains various options which are especially important to do nonlinear analysis. The nonlinear extensions to

466 MD Demonstration Problems CHAPTER 27

lower-order solid element, CHEXA can be activated by using the PSLDN1 property option to the regular PSOLID property option in the manner shown below: PSOLID PSLDN1 + C4

1 1

1 1 SOLI

0 1 L

+ +

The PLSLDN1 option allows the element to be used in both large displacement and large strain analysis and has no restrictions on the kinematics of deformation unlike the regular CHEXA elements with only PSOLID property entry. The standard CHEXA elements are more suitable for large rotations but small strain analysis due to their linear formulation in co-rotational system. While the difference may be small or even negligible in elastic analysis, use of any inelastic material model would certainly require the use of these options.

Modeling Contact The BCPARA defines the number of bodies in contact with maximum number of contact entities (e.g., patches), nodes on the periphery of the contact surfaces and contact parameters like friction type (in this case – node based, bilinear Coulomb model), friction coefficient, bias factor, and type of contact procedure used. BCPARA

0ERROR

0.005BIAS

0.99FTYPE

6

It must be mentioned that the contact procedure being used (flagged via ISPLIT flag) is iterative penetration checking procedure and must always be used for robustness in a quasi-static analysis. Friction has been flagged via the FTYPE field where a 6 denotes the bilinear, Coulomb model. The friction coefficient is 0.1 and is included in contact body definition with BCBODY option or the contact tables using the BCTABLE option. Another significant point is the use of BIAS in frictional problems. The bias factor measures the non-dimensionalized distance on both sides of the contact surface which is used to make a decision if the node is in contact or not, based on whether the node falls within this band defined by contact zone tolerance. Ideally, it should be 1.0 or as close to it. However, due to the possibility of excessive iterations in case of even very slight penetration, the bias is kept as zero or, in other words, a slight penetration is accepted. While a bias of zero works well for nonfrictional problems, it can be a detriment for frictional problems which require the bias to be set as close to one as possible in order to avoid a fictitious tangential force on the node which can cause non convergence of the solution. Finally, the ERROR parameter denotes the contact zone tolerance. The default value is about 1/20th of the smallest element size for a solid element. In this case, it has been chosen to be an even smaller value of 0.005. To identify how the contact bodies can touch each other, the BCTABLE option is used. BCTABLE with ID 0 is used to define the touching conditions at the start of the analysis. This is a mandatory option required in SOL 400 for contact analysis and it is flagged in the case control section through the optional BCONTACT = 0 option. The BCTABLE with ID 1 is used to define the touching conditions for later increments in the analysis, and it is flagged using BCONTACT = 1 in the case control section. Also, the SLAVE-MASTER combination defines that the nodes for body 1 are nodes belonging to the slave body. This, in literature, is referred by various terminologies as either contacting body nodes or tied nodes (imagining the situation of multi-point constraints). The nodes belonging to body 2 are said to belong to the master body which are also referred to as the contacted body nodes or the retained nodes (imagining the situation of multi-point constraints) BCTABLE

0

1 SLAVE MASTERS 1

2

0. 0

0. 0

.1 0

0.

0

CHAPTER 27 467 Large Sliding Contact Analysis of a Buckle

BCTABLE

1

1 SLAVE

2

0. 0

0. 0

.1 0

0.

0

MASTERS 1

The definition of the contact bodies (defined as body 1 and 2 in Figure 27-1) consists of the bulk data entries. The BCBODY option defines the deformable body including the body ID, dimensionality, type of body, type of contact constraints and friction etc. while the BSURF identifies the elements forming a part of the deformable body as: BCBODY BSURF

2 3D DEFORM 2 2 50000 50001 50002 50003 50007 50008 50009 50010 50011 50015 50016 50017 50018 50019 50023 50024 50025 50026 50027 … (list of element forming this body)

2 50004 50012 50020 50028

50005 50013 50021 50029

50006 50014 50022 50030

Material Modeling The isotropic, Hookean elastic material properties of the deformable body are defined using the following MAT1 option as follows: MAT1

1

10000.

0.4

Isotropi

The Young’s modulus is taken to be 10 GPa with a Poisson’s ratio of 0.4.

Loading and Boundary Conditions The displacements for body 2 are fixed at the end in the following manner: $ Displacement Constraints of Load Set : right_fixed_xyz SPC1 5 123 100056 THRU 100074 SPC1 5 123 100446 THRU 100464

The loading involves application of displacement controlled boundary conditions as follows: SPCADD 2 1 8 5 $ Enforced Displacements for Load Set : case1_left_xyz $ Dummy Force Required to Activate the Following Enforced Displacements FORCE 1 50084 0. .57735 .57735 .57735 SPCD 1 50084 1 85. 50085 1 85.

A total X displacement of 85 mm is applied to body 1. The application of the loads or displacements is such that the total load applied at the end of the loading sequence is given in the input.

Solution Procedure The nonlinear procedure used is defined through the following NLPARM entry: NLPARM NLAUTO

1 0.01 1 10

20 .01 0

FNT 1.

.1

1.2

50

UV

ALL

1.-5

.5

0

FNT represents Full Newton-Raphson technique wherein the stiffness is reformed at every iteration; KSTEP (field after FNT) is left blank and in conjunction with FNT, it indicates that the program will determine if the stiffness needs

468 MD Demonstration Problems CHAPTER 27

to be reformed between the end of the load step and the start of next load increment. Fifty (50) is the maximum number of allowed recycles for every increment and, if this were to be exceeded, the load step would be cut-back and the increment repeated. UV indicates that the maximum norm of vector component of the incremental displacements will be checked for convergence. ALL indicates that intermediate output will be produced after every increment. The second line of NLPARM indicates that a tolerance of 0.01 will be used for displacement based convergence checking. NLAUTO defines the parameters in the adaptive load stepping scheme. The initial load step is 1% of the total load. It must be noted that, for many problems including plasticity of complicated contact conditions in the early stages of the analysis, this must be a very small percentage (typically 0.5%). The smallest and largest ratio between the steps is 0.1 and 1.2, respectively, while the minimum value of the step is 10 – 5 . Finally, the desired number of recycles is kept at ten which is the default in SOL 400. If this number is chosen to be very small, then the step size is cut to a smaller size for convergence to be achieved and there will be larger number of steps. If this number is very large, then the load step will allow more iterations for convergence in the same step.

The number of increments is provided in the third field of the NLPARM option. It is also worth noting that removing the NLAUTO option results in a constant load step procedure with a total of 20 load increments per step (thus, a total of 40 for the analysis). Alternately another nonlinear procedure used is defined through the following NLSTEP entry like: NLSTEP + + +

1 ADAPT MECH

1. 1.00E-2 0 PV

1.E-5 0.0002

0.10 0.1

1.2

0

999999

+ + +

PFNT

Adaptive time procedure with total time of 1 is used. Initial time step of 0.01 is used as fraction of total time. It means the initial load step is 1% of the total load. It must be noted that, for many problems including plasticity of complicated contact conditions in the early stages of the analysis, this must be a very small percentage (typically 0.5%). The maximum number of recycles allowed for each increment are 10 and minimum is 1. The desired number of recycles per increment is 4. If this number is chosen to be very small, then the step size is cut to a smaller size for convergence to be achieved and there will be larger number of steps. If this number is very large, then the load step will allow more iterations for convergence in the same step.The smallest and largest ratio between the steps is 0.1 and 1.2, respectively, while the minimum value of the step is 1E-5. Output is written to result file for every single increment.

CHAPTER 27 469 Large Sliding Contact Analysis of a Buckle

Results Figure 27-2 shows the sequence of the analysis with a close-up view of the buckle. It can be seen that the clip slides on top of the protrusion of the static frame without any penetration. It is quite remarkable that even with the large motion as well as large sliding contact per load increment between the two deformable contact bodies, the analysis shows a robust behavior. A vector plot of the comparison of normal and frictional contact forces with the Marc results is presented in Figure 27-3 and Figure 27-4, respectively. The contact forces for SOL 400 and Marc agree very well in both magnitude and direction.

Figure 27-2

Various Stages of Insertion of the Clip

470 MD Demonstration Problems CHAPTER 27

(a) SOL 400

Figure 27-3

Comparison of Contact Normal Forces

(a) SOL 400

Figure 27-4

(b) Marc

(b) Marc

Comparison of Contact Frictional Forces

Next, the load displacement for the frictional and frictionless cases are compared in Figure 27-5. Only the X direction forces are plotted versus time. It is always recommended to perform a frictionless analysis (nug_27f.dat) whenever possible to aid in the understanding of the affect of adding friction. As expected, for the frictionless case, the load displacement curve is symmetric about the center line (between the insertion and removal steps). Deformed geometry is shown at various peaks of the curve and, as intuition would suggest, the peak forces correspond to the point of maximum bending. Addition of the non-conservative friction forces destroys the symmetry and the peak insertion force increases compared to the peak force in removal. The removal of the clip generates less pull-out force compared to the push-in force. Also, the insertion force starts reducing due to frictional forces aiding the motion as opposed to resisting the motion as the sliding switches from the convex part to the concave part of the contact surface.

CHAPTER 27 471 Large Sliding Contact Analysis of a Buckle

Fx

Fx

1000

Fx (N)

500 0.5

0

1.0

1.5 Time (s)

-500

Frictionless Frictional

-1000 -1500

Insert

Remove

-2000 Fx

Figure 27-5

Fx

Load Displacement Curve for the Frictional and Frictionless Cases

Checking the finite element analysis with a hand calculation assists both in understanding the FEM as well as the E t physics of the simulation. Solving elementary equations mentioned earlier for the bending stress yields,  = 3--- ----------2 2 L

where  is the tip displacement shown in Figure 27-6 during the insertion of the clip. Inc: 17 Time: 4.250e-001

4.213e+002

L=8

3.368e+002

0 mm

2.524e+002 1.679e+002 δ = 20 mm

8.349e+001 -9.664e-001

2ζ t = 6 mm

-8.542e+001 -1.699e+002 -2.543e+002 -3.388e+002 Y

-4.232e+002

lcase1 Comp 11 of Stress

Figure 27-6

Z

X 1

Verify FEM with Simple Calculation

Performing the calculation of the bending stress at the outer fibers of the thinnest section gives, 2 2 3  10x10 9  N  m    20mm   6mm  m N 3 E  2 t  N - =  4.69 x10 8 ------  ------------------ =  469 ----------- .  =  --- -------------------=  --- ----------------------------------------------------------------------------------2 2  10 3 mm 2 2 2 L2  80mm  m mm

The value of  469N   mm 2 

agrees closely to the corresponding bending stresses in Figure 27-6 of  423N   mm 2  . As expected, the linear solution presents an upper-bound to the actual stresses.

472 MD Demonstration Problems CHAPTER 27

Modeling Tips The two most important aspects in the analysis comprise of the inclusion of assumed strain enhancements to the standard element formulation and the choice of contact and time stepping scheme parameters use of adaptive load stepping scheme, and its associated parameters. The former is important due to presence of bending stresses in the structure which can manifest themselves as (sometimes large) spurious shear stresses. This is a purely numerical artifact due to the standard, displacement based finite element chosen which can be ameliorated by the use of an assumed strain enhancement to the standard element. Among the numerical parameters affecting the convergence of the job, the two most important parameters for this kind of analysis are the contact bias and maximum number of recycles for the adaptive stepping scheme. In contact analysis with friction, it is important to use a high bias (preferably 0.99) for frictional problems for improved convergent results. In many cases (although, not in this problem, nug_27b.dat), it can decrease the number of iterations as well. Next is the contact zone tolerance. Typically, a default value is 1/20th the smallest length of solid element. If the contact zone is too big, then there could be a loss of accuracy due to acceptance of penetrated nodes or large amount of recycling due to contact nodes separating. However, reducing the contact zone tolerance may not always yield the reduction in the number of iterations. In fact, in certain problems where there are not many separations expected, reducing to a very small number can even increase the number of iterations due to contact detection and scaling of incremental displacements in the iterative penetration checking algorithm in contact. It is also worth noting that the adaptive load stepping improves the speed and accuracy of the analysis quite significantly for this problem due to its intelligent choice of time steps based on the convergence parameters. This adequately demonstrates the strength of the adaptive stepping in tough problems where the smart algorithm adjusts the increment size based on the kinematics of deformation, contact constraints, and convergence rates rather than the fixed time stepping where the only alternative is to cut down the existing increment size in case of non convergence in the specified number of recycles. It is also noted that a very high or very low number of desired number of recycles can either invoke an excessive number of iterations or induce cutbacks during the analysis. For example, decreasing the desired number of recycles to may increase the number of increments. Due to a large amount of sliding and significant contact nonlinearity, a large number of recycles, in general, are expected for most increments. Therefore, a high number of desired recycles proved to be useful in this particular example. However, in problems with milder material and/or contact nonlinearities where only a few iterations per increment are expected, a smaller number of desired recycles can yield faster results. This difference can result in notable savings of the computing time for large jobs. Flat rigid surfaces can be glued to the ends of the buckle and insert to control the insertion and extraction of the insert in and out of the buckle. The advantage of this modeling technique is that the total insertion and extraction force component, Fx, can be easily determined as shown in Figure 27-5, since all of the forces acting on rigid bodies are resolved to a single force and moment vector acting at the position of the rigid bodies. Finally, since the buckle has a plane of symmetry, it is cost effective to only model the half of the model say above this plane of symmetry. Note:

For contact problems, artificial damping can improve the speed of convergence and stability of the analysis as seen in nug_27c.dat.

CHAPTER 27 473 Large Sliding Contact Analysis of a Buckle

Input File(s) File

Description

mug_27.dat

Marc input for fixed time

nug_27.dat

MD Nastran input for fixed time stepping

nug_27a.dat

MD Nastran input with adaptive time stepping with bias = 0.99, contact zone tolerance = 0.0 (default), desired number of recycles = 20 (default = 10)

nug_27b.dat

MD Nastran input with adaptive time stepping bias = 0.0 (default), contact zone tolerance = 0.005, desired number of recycles = 20 (default = 10)

nug_27c.dat

MD Nastran input with adaptive time stepping bias = 0.99, contact zone tolerance = 0.005, desired number of recycles = 20

nug_27b.bdf

Input file similar to nug_27b.dat above with half symmetry use in the video

nug_27_star t.SimXpert

MD Nastran input with adaptive time stepping bias = 0.99, contact zone tolerance = 0.005, desired number of recycles = 20

Video Click on the image or caption below to view a streaming video of this problem; it lasts approximately 47 minutes and explains how the steps are performed.

168 mm

Y

X Z

247

Figure 27-7

mm

y metr Sym Half

Video of the Above Steps