COMSOL 6.4 - Necking of an Elastoplastic Metal Bar

A circular bar of metal is subjected to a uniaxial tensile test. The material has an elastoplastic behavior with nonlinear isotropic hardening. When subjected to large deformations, the specimen experiences a significant plastic deformation and necking in its central cross section. This example demonstrates the large strain formulation of plasticity available in the Nonlinear Structural Materials Module. The simulation results are compared with results in the literature.

Model Definition

In this model, a cylindrical steel bar with a height, H0, of 53.334 mm and a radius, R0, of 6.413 mm, is subjected to a total elongation of 20 mm.

The problem exhibits 2D axial symmetry as well as a reflection symmetry in the mid cross section of the bar. It is therefore possible to reduce the model geometry to a rectangle with a width equal to the radius of the bar, and a height equal to half of the length of the bar, see Figure 1.

The boundary conditions for the displacements in the radial and axial directions, u and w, that follow from the symmetries are

where r and z are the radial and axial coordinates. The tensile loading is imposed through a prescribed elongation of Δ = 10 mm,

Figure 1: Schematic description of the numerical model. The red line shows the rotation symmetry axis.

Material Model

The elastic behavior of the material is characterized by a Young’s modulus E = 206.9 GPa and a Poisson’s ratio ν = 0.29. The plastic response follows a nonlinear isotropic hardening model with a yield stress given by

(1)

where σys0 is the initial yield stress and σh is the nonlinear hardening function. The latter is defined as

(2)

The hardening function σh depends nonlinearly on the equivalent plastic strain εpe. Here, H is the linear hardening coefficient, σysf is the saturation flow stress or residual yield stress, and ζ is the saturation exponent.

The numerical values of the parameters for the hardening function are given in Table 1. Figure 2, shows this nonlinear hardening as a function of equivalent plastic strain.

Table 1: Hardening function constants.
Constant	Value
σys0	450 MPa
H	129.24 MPa
σysf	715 MPa
ζ	16.93

Figure 2: Nonlinear isotropic hardening as a function of equivalent plastic strain.

As plastic strain become large, the material becomes unstable which will result in localization of the deformation in to the smallest possible width; in the case of finite elements, rows of mesh element or rows of Gauss points. To keep the mesh objectivity of the solution, a length scale can be introduced to the material model to limit strain localization to a predefined width. Here the nonlocal plasticity model is used with a length scale lint = 0.15 mm.

Results and Discussion

Figure 3 shows the distribution of the von Mises stress in the specimen. The figure also shows the deformation of the specimen, and how the nonlinear material behavior causes necking of its central cross section. The necking is caused by the very large plastic strains at the vertical symmetry plane as seen in Figure 4. A 10 mm elongation results in an equivalent plastic strain larger than 2.5 (this is, 250%) in the most affected region. At such large plastic strains, the specimen has likely already ruptured.

Figure 3: Distribution of the von Mises stress at 10 mm end displacement.

Figure 4: Distribution of the equivalent plastic strain at 10 mm end displacement. The results are shown in the undeformed configuration of the bar.

Figure 5 shows the change in radius as a function the elongation. Initially, the radius decreases linearly with the applied displacement. After an axial displacement of 3 mm, the radial reduction increases significantly and the specimen experiences necking. A similar case has been examined by Simo and Hughes (Ref. 1) as well as by Elguedj and Hughes (Ref. 2). Table 2 compares their results with the outcome of the current analysis. The calculated radial neck radius is in good agreement with results found in the literature.

Figure 5: Necking development in the mid section of the bar.

Table 2: Comparison between applied displacement and bar radius (mm) at the mid section.
δ (mm)	COMSOL	Ref. 1	Ref. 2
1.0	6.3	6.3	6.3
2.0	6.2	6.1	6.1
3.0	6.0	5.9	5.9
4.0	5.6	5.3	5.4
5.0	4.9	4.6	4.6
6.0	4.0	3.7	3.7

Figure 6 shows the reaction force as function of the axial displacement. The peak load is recorded at 3 mm after which the load starts to decrease as the specimen experience necking as discussed above. At an axial displacement around 7 mm, the curve starts to flatten, this is an effect of the nonlocal plasticity model. However, at this point the specimen in reality would have ruptured, but to capture such behavior, the material model would have to be extended.

Figure 6: Reaction force of the bar.

Notes About the COMSOL Implementation

The nonlinear hardening behavior is implemented using an analytic function. With reference to equation Equation 2, the hardening function is defined as

The user-defined hardening σh(εpe), portrayed in Figure 2, assumes zero stress for no plastic strain. In Equation 1, the yield stress is defined as the sum of the initial yield stress σys0 and the hardening function σh(εpe). The hardening function is defined in the node by an analytic function.

The large strain formulation of plasticity is used when the parent material model to he node requests a multiplicative decomposition of deformation gradients. Selecting the formulation in the section of will typically enable the multiplicative decomposition. By also setting to , one can force the multiplicative decomposition to be used always. As a rule-of-thumb, the multiplicative formulation is recommended when elastic and/or plastic strains exceed 0.1 (this is, more than 10%).

Including plasticity in the material model requires local computations at each Gauss point during the assembly process which is expensive. By using , the number of Gauss points are reduced by a factor two for the given displacement shape function order and mesh element type. This will speed up the computation significantly.

In this example, you use the double dogleg solver. This solver often works better for this class of nonlinear problems. The default Newton solver could also be used, but would then require tuning of the default solver settings.

References

1. J.C. Simo and T.R.J. Hughes, Computational Inelasticity, Springer, 2000.

2. T. Elguedj and T.J.R. Hughes, Isogeometric Analysis of Nearly Incompressible Large Strain Plasticity, ICES REPORT 11-35, The Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2011.

Nonlinear_Structural_Materials_Module/Plasticity/bar_necking

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
sigma0	450[MPa]	4.5E8 Pa	Initial yield stress
sigmaSF	715[MPa]	7.15E8 Pa	Saturation flow stress
H	129.24[MPa]	1.2924E8 Pa	Linear hardening coefficient
zeta	16.93	16.93	Saturation exponent
delta	0[m]	0 m	Top displacement
H0	53.334[mm]	0.053334 m	Bar length
R0	6.413[mm]	0.006413 m	Bar radius

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0.

In the text field, type H0/2.

Click

Solid Mechanics (solid)

Prescribed Displacement, Bottom

In the toolbar, click

and choose .

In the window for , type Prescribed Displacement, Bottom in the text field.

Select Boundary 2 only.

Locate the section. From the list, choose .

Prescribed Displacement, Top

In the toolbar, click

and choose .

In the window for , type Prescribed Displacement, Top in the text field.

Select Boundary 3 only.

Locate the section. From the list, choose .

From the list, choose .

In the u 0 z text field, type delta.

Linear Elastic Material 1

Use a large strain formulation with a multiplicative decomposition of deformation gradients.

In the window, click .

In the window for , locate the section.

From the list, choose .

The local computations at Gauss points during assembly are expensive when plasticity is added to the material model. Using a reduced integration scheme will reduce the overall simulation time.

Locate the section. Select the checkbox.

Plasticity 1

In the toolbar, click

and choose .

In the window for , locate the section.

Find the subsection. From the list, choose .

Click to expand the section. From the list, choose .

In the lint text field, type 0.15[mm].

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	206.9[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.29	1	Young’s modulus and Poisson’s ratio
Density	rho	7850	kg/m³	Basic
Initial yield stress	sigmags	sigma0	Pa	Elastoplastic material model