3.4 Example: stress wave propagation in a bar

This example demonstrates some of the fundamental ideas in explicit dynamics described earlier in the discussion of how stress propagates along the length of a bar. It also illustrates stability limits and the effect of mesh refinement and materials on the solution time.

The bar has the dimensions shown in Figure 3–5.

Figure 3–5 Problem description for wave propagation in a bar.

To make the problem a one-dimensional strain problem, all four lateral faces are on rollers; thus, the three-dimensional model simulates a one-dimensional problem. The material is steel with the properties shown in Figure 3–5. The free end of the bar is subjected to a blast load with a magnitude of 1.0 × 105 Pa, as shown in Figure 3–6, in which the duration of the blast load is 3.88 × 10–5 s.

Figure 3–6 Blast amplitude versus time.

Using the material properties (neglecting Poisson's ratio), we can calculate the wave speed of the material using the equations introduced in the previous section.

At this speed the wave passes to the fixed end of the bar in 1.94 × 10–4 s. Since we would like to see the stress propagate along the length of the bar through time, we need an adequately refined mesh to capture the stress wave accurately. It seems adequate to have the blast load take place over the span of 10 elements. Since the blast duration is 3.88 × 10–5 s, this means that we would like the blast duration times the wave speed to equal the length of 10 elements:

The length of 10 elements is 0.2 m. Since the length of the bar is 1.0 m, this means that we would like to have 50 elements along the length. To keep the mesh uniform, we will also have 10 elements in each of the transverse directions, making the mesh 50 × 10 × 10. This mesh is shown in Figure 3–7.

Figure 3–7 50 × 10 × 10 mesh.

Create this mesh in your preprocessor. Use the coordinate system shown in Figure 3–7.


3.4.1 Node and element sets

You will need to have node and element sets defined to apply the loads and boundary conditions and to visualize output. Define the node sets on their respective faces, as shown in Figure 3–8.

Figure 3–8 Node sets.

Define the element sets shown in Figure 3–9.

Figure 3–9 Element sets for modeling.

In addition, define an element set containing three elements in the center of the bar. Select these elements such that their faces nearest to the free end are at distances 0.25 m, 0.5 m, and 0.75 m from the free end, as shown in Figure 3–10. These elements will be used for postprocessing.

Figure 3–10 Element sets for postprocessing.


3.4.2 Reviewing the input file—the model data

At this point we assume that you have created the basic mesh using your preprocessor. In this section you will review your input file and include additional information.

Model description

The following would be a suitable description in the *HEADING option for this simulation:

*HEADING
Stress wave propagation in a bar -- 50x10x10 elements
SI units (kg, m, s, N)

Element connectivity

Check to make sure that you are using the correct element type (C3D8R). It is possible that the preprocessor specified the element type incorrectly. The *ELEMENT option block in your model should begin with the following:

*ELEMENT, TYPE=C3D8R, ELSET=BAR

The name given for the ELSET parameter in your model may not be BAR. If necessary, change the name to BAR.

Section properties

The section properties are the same for all elements. Use element set BAR to assign the material properties to the elements.

*SOLID SECTION, ELSET=BAR, MATERIAL=STEEL

Material properties

The bar is made of steel, which we assume to be linear elastic with a Young's modulus of 207 × 109 Pa, a Poisson's ratio of 0.3, and a density of 7800 kg/m3. Use the following material option block in your model:

*MATERIAL, NAME=STEEL
*ELASTIC
207.0E9, 0.3
*DENSITY
7800.0,

Fixed boundary conditions

Fix all the translations at the built-in, right-hand end of the bar. Constrain the front, back, top, and bottom faces of the bar so that these faces are on rollers and the strain is uniaxial. Using the node sets that you defined previously, apply the following boundary conditions to your model:

*BOUNDARY
NFIX, 1, 3
NFRONT, 3, 3
NBACK, 3, 3
NTOP, 2, 2
NBOT, 2, 2

Amplitude definition

The blast load is applied at its maximum value instantaneously and is held constant for 3.88 × 10–5 s. Then the load is suddenly removed and held constant at zero. The *AMPLITUDE option is used to define the time variation of loads and boundary conditions. On the data lines following the *AMPLITUDE option, pairs of data are given in the form:

<time>, <amplitude>, <time>, <amplitude>, etc.
Up to four data pairs can be entered on each data line. ABAQUS considers the amplitude to be held constant following the last amplitude value given. Define the amplitude for the blast load as follows:
*AMPLITUDE, NAME=BLAST
0., 1., 3.88E-5, 1., 3.89E-5, 0, 3.90E-5, 0.


3.4.3 Reviewing the input file—the history data

We will now review the history data associated with this problem, including the step definition, loading, bulk viscosity, and output requests.

Step definition

The step definition indicates that this is an explicit dynamics analysis with a duration of 2.0 × 10–4 s. You can also include a descriptive title for the step.

*STEP
Blast loading
*DYNAMIC, EXPLICIT
, 2.0E-4

Loading

Apply the pressure load with a value of 1.0 × 105 Pa to the free face of the bar, which you previously defined to be in an element set called ELOAD. The pressure load at any given time is the magnitude specified under the *DLOAD option times the value interpolated from the amplitude curve. To apply the load correctly, you need to determine the face identifier label of the free element faces. For the model defined in Stress wave propagation in a bar, Section A.2, the free face is face number 3, which corresponds to the pressure identifier P3. The face identifier depends on the order in which the nodes are defined on the *ELEMENT option, as shown in Figure 3–11. Use the amplitude named BLAST when applying the pressure load.

*DLOAD, AMPLITUDE=BLAST
ELOAD, <P1, P2, P3, P4, P5, or P6>, 1.0E5

If you define the pressure load in your preprocessor, the correct face identifier should be determined automatically.

Figure 3–11 Face label identifier for a C3D8R element.

Bulk viscosity

To keep the stress wave as sharp as possible, set the quadratic bulk viscosity (discussed in Bulk viscosity, Section 3.6.1) to zero.

*BULK VISCOSITY
, 0.0

Output requests

By default, many preprocessors create an ABAQUS input file that has a large number of output request options. If, when you edit your input file, you find that these default output options were created, delete them all because they will generally generate too much output.

You want to have an output database file created during the analysis so that you can use ABAQUS/Viewer to postprocess the results. Four output database frames (intervals at which data are written to the output database) are adequate to show the stress wave propagating through the mesh. Set the parameter VARIABLE=PRESELECT on the *OUTPUT, FIELD option to write the default field data for a *DYNAMIC, EXPLICIT procedure to the output database file. In addition, request stress (S) history output in element set EOUT every 1.0 × 10–6 s.

*OUTPUT, FIELD, VARIABLE=PRESELECT, NUMBER INTERVAL=4
*OUTPUT, HISTORY, TIME INTERVAL=1.0E-6
*ELEMENT OUTPUT, ELSET=EOUT 
S,
*END STEP


3.4.4 Running the analysis

After storing your input in a file called wave_50x10x10.inp, run the analysis using the following command:

abaqus job=wave_50x10x10
If your analysis does not complete, check the data file, wave_50x10x10.dat, and status file, wave_50x10x10.sta, for error messages. Modify your input file to remove the errors. If you still have trouble running your analysis, compare your input file to the one given in Stress wave propagation in a bar, Section A.2.

Status file

The status file, wave_50x10x10.sta, contains information about moments of inertia, followed by information concerning the initial stability limit. The 10 elements with the lowest stable time limits are listed in rank order.

    Most critical elements :
     Element number   Rank    Time increment   Increment ratio
    ----------------------------------------------------------
          101          1       1.931897E-06      1.000000E+00
          105          2       1.931897E-06      1.000000E+00
          107          3       1.931897E-06      1.000000E+00
          111          4       1.931897E-06      1.000000E+00
          115          5       1.931897E-06      1.000000E+00
          119          6       1.931897E-06      1.000000E+00
          401          7       1.931897E-06      1.000000E+00
          405          8       1.931897E-06      1.000000E+00
          407          9       1.931897E-06      1.000000E+00

The status file continues with information about the progress of the solution.

  STEP 1  ORIGIN 0.00000E+00
  
   Total memory used for step 1 is approximately 2.8 megawords
   Global time estimation algorithm will be used.
   Scaling factor :  1.0000
  
             STEP      TOTAL       CPU     STABLE      CRITICAL  KINETIC
  INCREMENT  TIME      TIME        TIME   INCREMENT    ELEMENT   ENERGY
         0  0.000E+00 0.000E+00   00:00:02 1.932E-06     101   0.000E+00
 Results number  0 at increment zero.
 ODB Field Frame Number    0 of    4 requested intervals at increment zero.
         6  1.159E-05 1.159E-05   00:00:05 1.932E-06     401   4.487E-05
        12  2.318E-05 2.318E-05   00:00:07 1.932E-06     101   9.566E-05
        16  3.332E-05 3.332E-05   00:00:09 3.104E-06     101   1.389E-04
        20  4.572E-05 4.572E-05   00:00:10 3.088E-06     101   1.712E-04
        22  5.189E-05 5.189E-05   00:00:11 3.082E-06     101   1.695E-04
 ODB Field Frame Number    1 of    4 requested intervals at  5.188920E-05
        26  6.420E-05 6.420E-05   00:00:13 3.071E-06     101   1.699E-04
.
.
.


3.4.5 Postprocessing

Start ABAQUS/Viewer by typing

abaqus viewer odb=wave_50x10x10
at the operating system prompt.

Plotting the stress along a path

We would like to look at how the stress distribution along the length of the bar changes with time. To do so, we will look at the stress distribution at three different times throughout the course of the analysis.

Create a curve of the stress in the 1-direction (S11) approximately along the center of the bar for each of the first three frames of the output database file. You first need to define a straight path along the center of the bar using a node list.

To create a node list path along the center of the bar:

  1. From the main menu bar, select ToolsPathCreate.

    The Create Path dialog box appears.

  2. Name the path Center. Accept the default selection of Node list as the path type, and click Continue.

    The Edit Node List Path dialog box appears.

  3. In the Node Labels table, enter the nodes along the center of the bar and click OK. For example, if you created your mesh using the input file shown in Stress wave propagation in a bar, Section A.2, the table entry is 100011:105011:100. (This input specifies a range of nodes from 100011 to 105011 in increments of 100.)

    Tip:  If you did not use the input file shown in Stress wave propagation in a bar, Section A.2, and need to identify the nodes along the center of the bar, plot the undeformed model shape and turn on the node labels using the undeformed plot options. From the main menu bar, select PlotUndeformed Shape; then select OptionsUndeformed Shape to open the Undeformed Shape Plot Options dialog box. You can use the view manipulation tools in the toolbar to alter the default isometric view of the bar.

To save X–Y plots of stress along the path at three different times:

  1. From the main menu bar, select ToolsXY DataManager.

  2. In the XY Data Manager dialog box that appears, click Create.

    The Create XY Data dialog box appears.

  3. Choose Path as the X–Y data source, and click Continue.

    The XY Data from Path dialog box appears with the path that you created visible in the list of available paths. If the undeformed model shape is currently displayed, the path you select is highlighted in the plot.

  4. Accept True distance as the selection in the X Values portion of the dialog box.

  5. Click Field Output in the Y Values portion of the dialog box to open the Field Output dialog box.

  6. Select the S11 stress component, and click OK.

    The field output variable in the XY Data from Path dialog box changes to indicate that stress data in the 1-direction (S11) will be created.

    Note:  ABAQUS/Viewer warns you that the field output variable will not affect the current image. Leave the plot mode As is, and click OK to continue.

  7. Click Step/Frame in the Y Values portion of the XY Data from Path dialog box.

  8. In the Step/Frame dialog box that appears, choose the second of the five recorded frames. (The first frame listed, frame 0, is the base state of the model at the beginning of the step.) Click OK.

    The Y Values portion of the XY Data from Path dialog box changes to indicate that data from Step 1, frame 1 will be created.

  9. To save the X–Y data, click Save As.

    The Save XY Data As dialog box appears.

  10. Name the X–Y data S11_T1, and click OK.

    S11_T1 appears in the XY Data Manager.

  11. Repeat Steps 7 through 10 to create X–Y data for frames 2 and 3. Name the data sets S11_T2 and S11_T3, respectively.

  12. To close the XY Data from Path dialog box, click Cancel.

To plot the stress curves:

  1. In the XY Data Manager dialog box, drag the cursor to highlight all three X–Y data sets.

  2. Click Plot.

    ABAQUS/Viewer plots the stress in the 1-direction along the center of the bar for frames 1, 2, and 3, corresponding to simulation times of 5.1208 × 10–5 s, 1.0016 × 10–4 s, and 1.5175 × 10–4 s, respectively.

To customize the X–Y plot:

  1. From the main menu bar, select OptionsXY Plot.

    The XY Plot Options dialog box appears.

  2. Click the Tick Marks tab.

    The Tick Marks options become available.

  3. Specify that the Y-axis major tick marks occur at 20 × 103 s increments.

  4. For both the X-axis and Y-axis minor tick mark options, specify zero minor tick marks between each major tick mark interval.

    You can also customize the axis titles.

  5. Click the Titles tab to make the Titles options available.

  6. In the X-axis area, select a User-specified title source. In the Title text field, enter Distance along bar.

  7. In the Y-axis area, specify Stress - S11 as the Y-axis title.

  8. Click OK to apply the customized X–Y plot parameters and to close the XY Plot Options dialog box.

To customize the appearance of the curves in the X–Y plot:

  1. From the main menu bar, select OptionsXY Curve.

    The XY Curve Options dialog box appears.

  2. In the XY Data field, select S11_T2.

  3. Choose the dotted line style for the S11_T2 curve, and click Apply.

    The S11_T2 curve becomes dotted.

  4. Repeat Steps 2 and 3 to make the S11_T3 curve dashed.

  5. Click Dismiss to close the XY Curve Options dialog box.

    The customized plot appears in Figure 3–12.

    Figure 3–12 Stress (S11) along the bar at three different times.

We can see that the length of the bar affected by the stress wave is approximately 0.2 m in each of the three curves. This distance should correspond to the distance that the blast wave travels during its time of application, which can be checked by a simple calculation. If the length of the wave front is 0.2 m and the wave speed is 5.15 × 103 m/s, the time it takes for the wave to travel 0.2 m is 3.88 × 10–5 s. As expected, this is the duration of the blast load that we applied. The stress wave is not exactly square as it passes along the bar. In particular, there is “ringing” or oscillation of the stress behind the sudden changes in stress. Linear bulk viscosity, discussed later in this chapter, damps the ringing so that it does not affect the results adversely.

Creating a history plot

Another way to study the results is to view the time history of stress at three different points within the bar.

To plot the stress history of the elements in set EOUT:

  1. From the main menu bar, select ResultHistory Output.

    ABAQUS/Viewer displays the History Output dialog box. The Output Variables field contains a list of all the variables in the history portion of the output database; these are the only variables you can plot. To see the complete description of the variable choices, increase the width of the History Output dialog box by dragging the right or left edge.

  2. Drag the cursor to highlight stress in the 1-direction (S11) for all three elements in element set EOUT.

  3. At the bottom of the History Output dialog box, click Plot.

    ABAQUS/Viewer displays an X–Y plot of the longitudinal stress in each element versus time.

  4. Click Dismiss to close the dialog box.

As before, you can customize the appearance of the plot.

To customize the X–Y plot:

  1. From the main menu bar, select OptionsXY Plot.

    The XY Plot Options dialog box appears.

  2. Click the Titles tab.

    The Titles options become available.

  3. In the X-axis area, specify Total time as the X-axis title.

  4. Click OK to apply the customized X–Y plot options and to close the dialog box.

To customize the appearance of the curves in the X–Y plot:

  1. From the main menu bar, select OptionsXY Curve.

    The XY Curve Options dialog box appears.

  2. In the XY Data field, select the temporary X–Y data label that corresponds to the element closest to the free end of the bar. (Of the elements in this set, this one is affected first by the stress wave.)

  3. Select a User-specified legend source.

  4. In the Legend text field, enter S11-0.25.

  5. Click Apply.

  6. In the XY Data field, select the temporary X–Y data label that corresponds to the element in the middle of the bar. (This is the element affected next by the stress wave.)

  7. Specify S11-0.5 as the curve legend text, and change the curve style to dotted.

  8. Click Apply.

  9. In the XY Data field, select the temporary X–Y data label that corresponds to the element closest to the fixed end of the bar. (This is the element affected last by the stress wave.)

  10. Specify S11-0.75 as the curve legend text, and change the curve style to dashed.

  11. Click OK to apply your changes and to close the dialog box.

    The customized plot appears in Figure 3–13.

    Figure 3–13 Time history of stress (S11) at three points along the length of the bar (0.25 m, 0.5 m, and 0.75 m).

In the history plot we can see that stress at a given point increases as the stress wave travels through the point. Once the stress wave has passed completely through the point, the stress at the point oscillates about zero.


3.4.6 How the mesh affects the stable time increment and CPU time

In Automatic time incrementation and stability, Section 3.3, we discussed how mesh refinement affects the stability limit and the CPU time. Here we will illustrate the effect with the wave propagation problem. We began with a reasonably refined mesh of square elements with 50 elements along the length and 10 elements in each of the two transverse directions. For illustrative purposes, we will now use a coarse mesh of 25 × 5 × 5 elements and observe how refining the mesh in the various directions changes the CPU time. The four meshes are shown in Figure 3–14.

Figure 3–14 Meshes from least to most refined.

Table 3–1 shows how the CPU time (normalized with respect to the coarse mesh model result) changes with mesh refinement for this problem. The first half of the table provides the expected results, based on the simplified stability equations presented in this guide; the second half of the table provides the results obtained by running the analyses in ABAQUS/Explicit on a desktop workstation.

Table 3–1 Mesh refinement and solution time.

MeshSimplified TheoryActual
(s)Number of ElementsCPU Time (s)Max (s)Number of ElementsNormalized CPU Time
25 × 5 × 5ABC6.06e-66251
50 × 5 × 5A/22B4C3.14e-612504
50 × 10 × 5A/24B8C3.12e-625008.33
50 × 10 × 10A/28B16C3.11e-6500016.67

For the theoretical results we choose the coarsest mesh, 25 × 5 × 5, as the base state, and we define the stable time increment, the number of elements, and the CPU time as variables A, B, and C, respectively. As the mesh is refined, two things happen: the smallest element dimension decreases, and the number of elements in the mesh increases. Each of these effects increases the CPU time. In the first level of refinement, the 50 × 5 × 5 mesh, the smallest element dimension is cut in half and the number of elements is doubled, increasing the CPU time by a factor of four over the previous mesh. However, further doubling the mesh to 50 × 10 × 5 does not change the smallest element dimension; it only doubles the number of elements. Therefore, the CPU time increases by only a factor of two over the 50 × 5 × 5 mesh. Further refining the mesh so that the elements are uniform and square in the 50 × 10 × 10 mesh again doubles the number of elements and the CPU time.

This simplified calculation predicts quite well the trends of how mesh refinement affects the stable time increment and CPU time. However, there are reasons why we did not compare the predicted and actual stable time increment values. First, recall that we made the approximation that the stable time increment is

We then assumed that the characteristic element length, , is the smallest element dimension, whereas ABAQUS/Explicit actually determines the characteristic element length based on the overall size and shape of the element. Another complication is that ABAQUS/Explicit employs a global stability estimator, which allows a larger stable time increment to be used. These factors make it difficult to predict the stable time increment accurately before running the analysis. However, since the trends follow nicely from the simplified theory, it is straightforward to predict how the stable time increment will change with mesh refinement.


3.4.7 How the material affects the stable time increment and CPU time

The same wave propagation analysis performed on different materials would take different amounts of computer time, depending on the wave speed of the material. For example, if we were to change the material from steel to aluminum, the wave speed would change from 5.15 × 103 m/s to

The change from aluminum to steel has minimal effect on the stable time increment, because the stiffness and the density differ by nearly the same amount. In the case of lead the difference is more substantial, as the wave speed decreases to

which is approximately one-fifth the wave speed of steel. The stable time increment for the lead bar would be five times the stable time increment of our steel bar.