Viscoelastic rod subjected to constant axial load

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

This page discusses:

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

ER(t)=k1+k2e-t/τR,

where τR=η/k2, η is the damping coefficient and k1 and k2 are constants. In this case k1 is 6.89 MPa (1000 psi); k2 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, E0, 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.

E0 is immediately available as E0=k1+k2= 68.9 MPa (10000 psi), and ν0 is

ν0=3K-k1-k26K=0.4833.

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

GR(t)=G0(1-g¯1P(1-e-t/τ1G)).

We need to compute G0, g¯1P, and τ1G 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) properties of the material approach that of a linear elastic solid, with E=k1. The long-term shear modulus, G, can be calculated using the relationship between the bulk, shear, and extensional moduli:

G=3Kk19K-k1.

Likewise, the “instantaneous” or “glassy” shear modulus, G0, is

G0=3K(k1+k2)9K-(k1+k2).

Then g¯1P=1-G/G0= 0.901001.

The shear relaxation time, τ1G, 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:

dGRdt|t=0=(dGRdERdERdt)|t=0.

After some algebra we obtain

τ1G=9K-k1-k29K-k1τR=0.9899sec.

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 C10= 6.89 MPa (1000 psi), C01= 4.59 MPa (666.67 psi), and D1= 1.378 × 10−7 MPa−1 (0.00002 psi−1). These constants are such that the initial Young's modulus and initial Poisson's ratio are equal to E0 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−3. 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:

u2L(t=0)=σ22k1+k2=0.01u2L(t)=σ22k1=0.1.

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:

ν(t=0)=3K-k1-k26K=0.4833ν(t)=3K-k16K=0.4983.

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:

ν(t=0.001)=-ε11ε22=--4.8353×10-31.0005×10-2=0.4833ν(t=50)=-ε11ε22=--4.9478×10-29.9291×10-2=0.4983.

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, ER(t), and the time-dependent creep compliance, JC(t):

0tJC(τ)ER(t-τ)dτ=t.

Differentiation of this relation with respect to t yields

JC(t)    ER(0)+0tJC(τ)E˙R(t-τ)dτ=1.

With the previously used expression for ER(t) this takes the form

(k1+k2)JC(t)-k2τR0tJC(τ)e(τ-t)/τRdτ=1.

Differentiating this expression once more provides

(k1+k2)J˙C(t)-k2τRJC(t)+k2τR20tJC(τ)e(τ-t)/τRdτ=0.

Multiplying this equation by τR and adding it to the previous equation yields the differential equation

(k1+k2)τRJ˙C(t)+k1JC(t)=1.

With the introduction of the creep time constant τC=(1+k2/k1)τR this can be written as:

τCJ˙C(t)+JC(t)=1/k1.

The general solution to this differential equation is

JC(t)=1/k1+Ce-t/τC,

where the coefficient C is defined by the initial condition JC(0)=1/(k1+k2), which yields

JC(t)=(1-k2k1+k2e-t/τC)1k1.

For this problem the stress σ22 is a constant, so that

ε22=JC(t)σ22=(1-k2k1+k2e-t/τC)σ22k1.

From the values given above for k1, k2, and τR, as well as the fact that σ22=0.689 MPa (100 psi), ε22 becomes

ε22=0.1(1-0.9e-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

  1. Collingwood G. A.EBBecker, and TMiller, User's Manual for the TEXVISC Computer Program, Morton Thiokol, Inc., Document Numbers U-85-4550A and U-85-4550B, 1985.
  2. 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.