Viscoelastic rod subjected to constant axial load

This example verifies the use of the time domain linear viscoelastic material model in Abaqus.

The following topics are discussed:

Problem description
Material
Loading
Results and discussion
Input files
References
Figures

ProductsAbaqus/Standard

Problem description

The example, taken from Collingwood et al. (1985), is a rod of length of 254 mm (10 in) and diameter of 25.4 mm (1 in). The rod is fixed in the axial direction on one end and a constant axial load of 0.689 MPa (100 psi) is applied suddenly to the other end. The rod is modeled using one quadratic, axisymmetric, hybrid continuum element (CAX8H).

Material

The linear viscoelastic material model used in this example can be represented by a combination of linear springs and a dashpot, as shown in Figure 1. The extensional relaxation function is

E_{R} (t) = k_{1} + k_{2} e^{- t / τ_{R}},

where $τ_{R} = η / k_{2}$ , $η$ is the damping coefficient and $k_{1}$ and $k_{2}$ are constants. In this case $k_{1}$ is 6.89 MPa (1000 psi); $k_{2}$ is 62.01 MPa (9000 psi); $τ_{R}$ is 1.0 sec; and the bulk modulus, $K$ , is 689 MPa (100,000 psi) and is independent of time.

Short-term material properties are specified using linear elastic behavior (Linear elastic behavior), which requires the instantaneous Young's modulus, $E_{0}$ , and Poisson's ratio, $ν_{0}$ . The time-dependent behavior is specified using viscoelastic behavior, in which the shear relaxation modulus and the bulk modulus are defined by a Prony series (see Time domain viscoelasticity). For the Abaqus analysis of this problem, it is assumed that no volumetric relaxation occurs.

$E_{0}$ is immediately available as $E_{0} = k_{1} + k_{2} =$ 68.9 MPa (10000 psi), and $ν_{0}$ is

ν_{0} = \frac{3 K - k_{1} - k_{2}}{6 K} = 0.4833.

The time-dependent material behavior is approximated with a single-term Prony series for the shear relaxation modulus:

G_{R} (t) = G_{0} (1 - {\bar{g}}_{1}^{P} (1 - e^{- t / τ_{1}^{G}})) .

We need to compute $G_{0}$ , ${\bar{g}}_{1}^{P},$ and $τ_{1}^{G}$ from the extensional relaxation function. The limiting cases of both the shear and extensional relaxation functions are used for this purpose. The long-term ( $t \to \infty$ ) properties of the material approach that of a linear elastic solid, with $E = k_{1}$ . The long-term shear modulus, G, can be calculated using the relationship between the bulk, shear, and extensional moduli:

G = \frac{3 K k_{1}}{9 K - k_{1}} .

Likewise, the “instantaneous” or “glassy” shear modulus, $G_{0}$ , is

G_{0} = \frac{3 K (k_{1} + k_{2})}{9 K - (k_{1} + k_{2})} .

Then ${\bar{g}}_{1}^{P} = 1 - G / G_{0} =$ 0.901001.

The shear relaxation time, $τ_{1}^{G}$ , is obtained by writing the rate of change of the shear modulus in terms of the rate of change of the extensional modulus at time $t =$ 0:

{\frac{d G_{R}}{d t} |}_{t = 0} = {(\frac{d G_{R}}{d E_{R}} \frac{d E_{R}}{d t}) |}_{t = 0} .

After some algebra we obtain

τ_{1}^{G} = \frac{9 K - k_{1} - k_{2}}{9 K - k_{1}} τ_{R} = 0.9899 sec .

The same problem is also treated as a large-strain example. The relaxation behavior is defined in the same way, but the short-term elastic properties are given using hyperelastic behavior. The polynomial formulation with $N =$ 1 is used, and the constants are $C_{10} =$ 6.89 MPa (1000 psi), $C_{01} =$ 4.59 MPa (666.67 psi), and $D_{1} =$ 1.378 × 10⁻⁷ MPa⁻¹ (0.00002 psi⁻¹). These constants are such that the initial Young's modulus and initial Poisson's ratio are equal to $E_{0}$ and $ν_{0}$ , respectively, and produce a close fit to a linear material. (See Hyperelastic behavior of rubberlike materials for further discussion of the choice of constants when $N =$ 1.)

Loading

A distributed load of 0.689 MPa (100 psi) is applied instantaneously and held constant throughout the analysis. To model this, we use the quasi-static procedure (Quasi-static analysis) in two steps. The load is applied in the first step, which has a time period of 0.001 seconds, so that the instantaneous (glassy) behavior dominates. Since this step uses only one increment, a tolerance to control the accuracy of the creep integration is not specified. The second step has a time period of 50 seconds, during which the load is held constant and the rod is allowed to relax toward its long-term behavior. Automatic time incrementation is chosen by giving a tolerance value for the maximum difference in the creep strain increment over a time increment. This tolerance is selected so that its value is of the same order of magnitude as the maximum elastic strain. Therefore, for this example the tolerance is set to 5 × 10⁻³. The total force and the total moment on the loaded face of the model are output to the results file.

Results and discussion

The instantaneous and long-term behaviors provide a check on the Abaqus results. The instantaneous and long-term axial displacements of the rod tip can be calculated as follows:

\begin{matrix} \frac{u_{2}}{L} (t = 0) & = \frac{σ_{22}}{k_{1} + k_{2}} = 0.01 \\ \frac{u_{2}}{L} (t \to \infty) & = \frac{σ_{22}}{k_{1}} = 0.1. \end{matrix}

These values agree well with the Abaqus results. Similarly, the instantaneous and long-term values of the Poisson's ratio can also be calculated exactly:

\begin{matrix} ν (t = 0) & = \frac{3 K - k_{1} - k_{2}}{6 K} = 0.4833 \\ ν (t \to \infty) & = \frac{3 K - k_{1}}{6 K} = 0.4983. \end{matrix}

The Poisson's ratio can be extracted from the Abaqus results by taking the ratio of the lateral strain to the axial strain at $t =$ 0.001 and $t =$ 50:

\begin{matrix} ν (t = 0.001) & = - \frac{ε_{11}}{ε_{22}} = - \frac{- 4.8353 \times 10^{- 3}}{1.0005 \times 10^{- 2}} = 0.4833 \\ ν (t = 50) & = - \frac{ε_{11}}{ε_{22}} = - \frac{- 4.9478 \times 10^{- 2}}{9.9291 \times 10^{- 2}} = 0.4983. \end{matrix}

Since this is an applied stress problem, obtaining the exact solution for the entire time period of the analysis requires inverting the original constitutive integral equation defining uniaxial stress in terms of uniaxial strain. To perform this inversion, we use the following relation (Pipkin, 1972) between the time-dependent relaxation modulus, $E_{R} (t)$ , and the time-dependent creep compliance, $J_{C} (t)$ :

\int_{0}^{t} J_{C} (τ) E_{R} (t - τ) d τ = t .

Differentiation of this relation with respect to t yields

J_{C} (t) E_{R} (0) + \int_{0}^{t} J_{C} (τ) {\dot{E}}_{R} (t - τ) d τ = 1.

With the previously used expression for $E_{R} (t)$ this takes the form

(k_{1} + k_{2}) J_{C} (t) - \frac{k_{2}}{τ_{R}} \int_{0}^{t} J_{C} (τ) e^{(τ - t) / τ_{R}} d τ = 1.

Differentiating this expression once more provides

(k_{1} + k_{2}) {\dot{J}}_{C} (t) - \frac{k_{2}}{τ_{R}} J_{C} (t) + \frac{k_{2}}{τ_{R}^{2}} \int_{0}^{t} J_{C} (τ) e^{(τ - t) / τ_{R}} d τ = 0.

Multiplying this equation by $τ_{R}$ and adding it to the previous equation yields the differential equation

(k_{1} + k_{2}) τ_{R} {\dot{J}}_{C} (t) + k_{1} J_{C} (t) = 1.

With the introduction of the creep time constant $τ_{C} = (1 + k_{2} / k_{1}) τ_{R}$ this can be written as:

τ_{C} {\dot{J}}_{C} (t) + J_{C} (t) = 1 / k_{1} .

The general solution to this differential equation is

J_{C} (t) = 1 / k_{1} + C e^{- t / τ_{C}},

where the coefficient C is defined by the initial condition $J_{C} (0) = 1 / (k_{1} + k_{2})$ , which yields

J_{C} (t) = (1 - \frac{k_{2}}{k_{1} + k_{2}} e^{- t / τ_{C}}) \frac{1}{k_{1}} .

For this problem the stress $σ_{22}$ is a constant, so that

ε_{22} = J_{C} (t) σ_{22} = (1 - \frac{k_{2}}{k_{1} + k_{2}} e^{- t / τ_{C}}) \frac{σ_{22}}{k_{1}} .

From the values given above for $k_{1}$ , $k_{2},$ and $τ_{R}$ , as well as the fact that $σ_{22} =$ 0.689 MPa (100 psi), $ε_{22}$ becomes

ε_{22} = 0.1 (1 - 0.9 e^{- t / 10}) .

From this equation, it is evident that the effective time constant for the problem is dramatically different (by a factor of 10 in this case), depending upon whether the loading is applied force or applied displacement. Figure 2 is a time history plot of $ε_{22}$ as predicted by the above equation and as calculated by Abaqus. The plot shows acceptable agreement between the Abaqus results and the exact solution. Closer agreement can be obtained by using a smaller tolerance value for the maximum difference in the creep strain increment over a time increment.

The solution obtained with the large-strain formulation differs negligibly from the small-strain solution. Abaqus automatically converts frequency domain data into a time domain Prony series representation. The analysis results using Prony parameters calibrated from tabulated frequency-dependent moduli data are in good agreement with the analyses using time domain data directly.

Input files

viscorod_smallstrain.inp: Small-strain input data for this problem.
viscorod_small_frq2tim.inp: Small-strain input data for time domain analysis using Prony parameters calibrated from tabulated frequency-dependent moduli data.
viscorod_largestrain.inp: Large-strain input data for this problem.
viscorod_large_frq2tim.inp: Large-strain input data for this problem using Prony parameters calibrated from tabulated frequency-dependent moduli data.
viscorod_c3d8.inp: Model using the three-dimensional 8-node brick element, C3D8.
viscorod_cps4.inp: Model using the 4-node plane stress element, CPS4.
viscorod_t3d2.inp: Model using the 2-node truss element, T3D2.
viscorod_postoutput.inp: POST OUTPUT job for the restart file generated in viscorod_largestrain.inp.

References

Collingwood, G. A., E. B. Becker, and T. Miller, User's Manual for the TEXVISC Computer Program, Morton Thiokol, Inc., Document Numbers U-85-4550A and U-85-4550B, 1985.
Pipkin, A. C., Lectures on Viscoelasticity Theory, Springer Verlag, New York, 1972.

Figures

Figure 1. Spring and dashpot model of viscoelastic rod.

Figure 2. Time history of strain in the direction of load for viscoelastic rod.