We will represent a cylindrical domain using our trapezoidal mesh. The axis of rotation is the left side of the trapezoid. We also use axial symmetry on the bottom surface to represent a roughly hourglass shaped domain that is being pulled on the bottom and the top.

The first change to make is the names of the displacements. Throughout the input file, change disp_x to disp_r and disp_y to disp_z. Now, add a new block defining the axisymmetric problem

  coord_type = RZ

Since we are using the tensor mechanics master action, we do not need to make any changes to the strain calculator; the tensor mechanics action automatically detects the coordinate system from coord_type and sets the appropriate strain calculator

Now that the model is axisymmetric, we are calculating the stress in the $$$\theta$$$ direction or the hoop stress. The hoop stress is the (2,2) component. To visualize this stress component, change the option for the stress tensor component in the master action from stress_xx to stress_zz. The master action uses a generic XYZ coordinate system to determine the output variables, regardless of the names of the displacement variables given.

You can now check your input file against part_2.1.i. When you run the input file, the hoop stress should look like this:

Hoop stress for 2D axisymmetry and linear elasticity

Finite Strain Simulation

Up until now, all of our simulations have assumed small and total strains, and we applied the strain all at once. Now, we will switch to a finite strain elastic formulation where we apply the strain in increments over time.

The first change we make is to the strain parameter in the tensor mechanics master action to use the FINITE option:

    strain = FINITE

Then we will make a corresponding change to the stress material block to use a stress calculator that operates on finite strain increments:

  type = ComputeFiniteStrainElasticStress

Because we will now apply the strain over increments rather than all at once, we also need to change our applied deformation. We will apply a deformation of $$$0.0007\times t$$$, where the time $$$t\leq5$$$. The fixed boundary conditions stay the same, but the top boundary condition becomes

  type = FunctionPresetBC
  variable = disp_z
  boundary = top
  function = '0.0007*t'

Finally, because we are now marching through time, we need to change from a steady state simulation to a transient one. Our simulation will end at $$$t = 5$$$ and will take time steps of 1 (which is the default). Also, we will switch to a preconditioned JFNK solve. Our Executioner block becomes

  type = Transient
  end_time = 5

  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap -ksp_gmres_restart'
  petsc_options_value = 'asm      lu           1               101'

To improve our accuracy we will switch to a second order Lagrange element. To do this, we need to make a change in the Mesh block. To change the mesh, add to the Mesh block

second_order = true

The tensor mechanics master action will automatically detect the second order mesh file and will set the displacement variables to use the desired second order Lagrange shape functions.

You can check your input file against part_2.2.i. Now, when you run the simulation, you will see the solve information for five time steps. The final deformation and stress should still look about the same as last time.

Larger Deformations with J2 Plasticity

Now we will change from only considering elastic deformation to including plastic deformation. We will represent the plastic deformation using J2 plasticity with a yield stress of 240. We will initially assume no hardening, i.e. the model will be elastic/perfectly plastic. The plastic deformation is included by changing our stress material block to ComputeMultiPlasticityStress. This material can consider multiple yield surfaces, however we will consider only one. The yield behavior is defined in a new block, the UserObject block. UserObjects are flexible tools for computing necessary information required by other objects. In this case, they are needed by the ComputeMultiPlasticityStress material. Change the stress block to this:

  type = ComputeMultiPlasticityStress
  ep_plastic_tolerance = 1e-9
  plastic_models = J2

Now, add the UserObjects block

    type = TensorMechanicsHardeningConstant
    value = 3.0e2
    type = TensorMechanicsPlasticJ2
    yield_strength = str
    yield_function_tolerance = 1E-3
    internal_constraint_tolerance = 1E-9

You can learn more about UserObjects by looking here.

Now, increase the end time from 5 to 20. Also, a smaller time step is needed to fully capture the plastic behavior, so make dt = 0.25. This simulation can now take a while to complete, so coarsen the mesh back to the original size by changing uniform_refine = 0 in the mesh block.

We also want to be able to plot a stress-strain curve from our simulation. To accomplish this, we need to output the (1,1) component, along the axisymetric z-axis, of both the stress and the strain. As before, we'll use the generate_output option in the tensor mechanics master action to generate these quantities. Recall that master action uses generic XYZ notation for tensor components, so we will use the _yy options for the (1,1) components of these two tensors.

Change the tensor mechanics master action generate_output parameter to:

    generate_output = 'stress_yy strain_yy'

To get a single value for the stress and the strain, we will average them over the domain using a postprocessor. Postprocessors calculate scalar values that provide information about the solve or the variables. In this case, we will use a postprocessor that calculates the average value along a boundary. Add this block to your input file:

    type = SideAverageValue
    variable = stress_yy
    boundary = bottom
    type = SideAverageValue
    variable = strain_yy
    boundary = bottom

Note that the variable names we use in the Postprocessors block match the options for the generate_output parameter in the tensor mechanics master action. To learn more about Postprocessors, look here.

Now run your simulations. You will see that the solves are taking many more linear and nonlinear iterations. To simplify the output, you can suppress the output of the linear residuals by adding to your Outputs block the line

print_linear_residuals = false

To plot the stress strain curve, we need to output the Postprocessor values to a csv file by adding to the Outputs block the line

csv = true

To check your input file, compare against part_2.3.i. Now, run the input file and plot a stress-strain curve from the data in the .csv file. The stress strain curve has a linear region, and then yields and becomes flat because of our assumptions of elastic-perfectly plastic behavior, as shown to the right.

Stress-strain curve for elastic/perfectly plastic behavior

Adding A Hardening Law

Real materials are typically not perfectly plastic. Once plastic deformation begins, hardening causes the yield stress to increase. We will using a cubic hardening model with an initial yield stress of 240 and a maximum yield stress of 300. To do this, we change the hardening law in the UserObjects block. Change the hardening block to

  type = TensorMechanicsHardeningCubic
  value_0 = 2.4e2
  value_residual = 3.0e2
  internal_0 = 0
  internal_limit = 0.005

Look at part_2.4.i to check your input file. Now run your simulation and again plot the stress strain curve. You can see the impact of hardening on the behavior of the uniaxial tension sample, as shown to the right.

Stress-strain curve for a 2D axisymmetric plastic problem with cubic hardening

Part 3: 3D Simulation