Time-dependent three-dimensional quasi-static analysis of a viscoelastic solid by defining a time function

Abstract

In this paper, we present a novel formulation for the time-dependent three-dimensional bending analysis of a viscoelastic solid by defining a time function with unknown coefficient. The advantage of this method is that there is no need to follow the time-displacement curve; instead, the time-dependent responses of the viscoelastic solid are obtained only by bending analysis of the elastic solid with low computational cost and high computational speed. The relationship between stress and strain is written for linear viscoelastic materials employing the Boltzmann integral law. The displacement field is approximated through the combination of two functions, a function of geometrical parameters and a function of time. The equilibrium discretized equations are written based on the virtual work principle using the finite element method. The numerical results are compared with other available references. To determine the effects of geometry and material on the coefficient of the time function, the results are extracted for square, rectangular, skew, circular, and triangular prismatic solids as well as for different materials.

Appendices

Appendix A

In this section, we formulate the quasi-static analysis of the circular viscoelastic prism. We use the axial symmetry for the 2D analysis of the solid with a circular cross-section.

1.1 A.1 Kinematics

As shown in Fig. 37, the 2D time-dependent displacement field of the circular prismatic viscoelastic solid in the cylindrical coordinate system is written as follows:

$$ \mathbf{U} \left ( r, z, t \right ) = \left [ \textstyle\begin{array}{c} u \left ( r, z, t \right )\\ w \left ( r, z, t \right ) \end{array}\displaystyle \right ], $$

(43)

where $u \left ( r, z, t \right )$ and $w \left ( r, z, t \right )$ are the displacements in the $r$ and $z$ directions, respectively.

Fig. 37

The circular viscoelastic prism

The 2D strain vector of the solid is expressed as follows:

$$ \boldsymbol{\varepsilon } \left ( r, z, t \right ) = \left [ \textstyle\begin{array}{c} \varepsilon _{r}\\ \varepsilon _{\theta } \\ \textstyle\begin{array}{c} \varepsilon _{z}\\ \gamma _{rz} \end{array}\displaystyle \end{array}\displaystyle \right ] = \left [ \textstyle\begin{array}{c} \textstyle\begin{array}{c@{\quad }c} \frac{\partial }{\partial r} & 0 \end{array}\displaystyle \\ \textstyle\begin{array}{c@{\quad }c} \frac{1}{r} & 0 \end{array}\displaystyle \\ \textstyle\begin{array}{c} \textstyle\begin{array}{c@{\quad }c} 0 & \frac{\partial }{\partial z} \end{array}\displaystyle \\ \textstyle\begin{array}{c@{\quad }c} \frac{\partial }{\partial z} & \frac{\partial }{\partial r} \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \right ] \left [ \textstyle\begin{array}{c} u \left ( r, z, t \right )\\ w \left ( r, z, t \right ) \end{array}\displaystyle \right ]. $$

(44)

1.2 A.2 The stress–strain relationship

We write the relaxation modulus tensor of the circular viscoelastic prism as follows:

$$ \mathbf{D} \left ( t \right ) = \frac{E \left ( t \right )}{\left ( 1+ \nu \left ( t \right ) \right ) \left ( 1-2 \nu \left ( t \right ) \right )} \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \textstyle\begin{array}{c} 1- \nu \left ( t \right )\\ \nu \left ( t \right )\\ \textstyle\begin{array}{c} \nu \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} \nu \left ( t \right )\\ 1- \nu \left ( t \right )\\ \textstyle\begin{array}{c} \nu \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c@{\quad }c} \textstyle\begin{array}{c} \nu \left ( t \right )\\ \nu \left ( t \right )\\ \textstyle\begin{array}{c} 1- \nu \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} 0\\ 0\\ \textstyle\begin{array}{c} 0\\ 0.5- \nu \left ( t \right ) \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \right ]. $$

(45)

By inserting Eqs. (12) into Eq. (45) we obtain the relaxation modulus tensor in terms of the relaxation function:

$$ \mathbf{D} \left ( t \right ) = K \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \textstyle\begin{array}{c} 1+2 \omega \left ( t \right )\\ 1- \omega \left ( t \right )\\ \textstyle\begin{array}{c} 1- \omega \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} 1- \omega \left ( t \right )\\ 1+2 \omega \left ( t \right )\\ \textstyle\begin{array}{c} 1- \omega \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c@{\quad }c} \textstyle\begin{array}{c} 1- \omega \left ( t \right )\\ 1- \omega \left ( t \right )\\ \textstyle\begin{array}{c} 1+2 \omega \left ( t \right )\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} 0\\ 0\\ \textstyle\begin{array}{c} 0\\ 1.5 \omega \left ( t \right ) \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \right ], $$

(46)

where $K$ is the constant bulk modulus.

1.3 A.3 Approximating the displacement field

As discussed in Sect. 2.3, we can approximate the displacement field using Eqs. (47) in the cylindrical coordinate system:

$$ \mathbf{U} \left ( r, z, t \right ) = \left [ \textstyle\begin{array}{c} u \left ( r, z, t \right )\\ w \left ( r, z, t \right ) \end{array}\displaystyle \right ] = \mathbf{U}^{0} \left ( r, z \right ) + \mathbf{U}^{c} \left ( r, z \right ) f \left ( t \right ),\qquad \mathbf{U}^{c} = \mathbf{U}^{\infty } - \mathbf{U}^{0}, $$

(47)

where $\mathbf{U}^{0}$ and $\mathbf{U}^{\infty }$ are the displacement vectors at $t =0$ and $t =\infty$, respectively, obtained through the bending analysis of the elastic solid using the material properties at $t =0$ ($E^{0} = E \left ( 0 \right )$, $\nu ^{0} = \nu \left ( 0 \right )$) and $t =\infty$ ($E^{\infty } = E \left ( \infty \right )$, $\nu ^{\infty } = \nu \left ( \infty \right )$).

1.4 A.4 Virtual work principle

For a solid subjected to the distributed load $q$, according to the virtual work principle, the equilibrium equation can be written as follows:

$$ {2 \pi } \int _{{A}} {\delta } \boldsymbol{\upvarepsilon }^{\boldsymbol{T}} \boldsymbol{\upsigma } {rdA} ={2 \pi } \int _{{R}} {{q\delta wrdr}}. $$

(48)

1.5 A.5 Discretizing equations

Based on the finite element relations, we write the displacement function of each element as follows:

$$ \mathbf{U}_{i} \left ( r, z, t \right ) = \mathbf{n} \left ( r, z \right ) \mathbf{u}_{i}^{0} + \mathbf{n} \left ( r, z \right ) \left ( \mathbf{u}_{i}^{\infty } - \mathbf{u}_{i}^{0} \right ) f \left ( t \right ), $$

(49)

where $\mathbf{u}_{i}^{\infty }$ and $\boldsymbol{u}_{i}^{0}$ are the displacement matrices of the $i$th element nodes at $t =\infty$ and $t =0$, respectively, and $\mathbf{n} \left ( r, z \right )$ represents the shape functions of the elements.

We write the displacement function of each 1D element as follows:

$$ w_{j} \left ( r, z, t \right ) = \mathbf{n}_{s} \left ( r, z \right ) \mathbf{w}_{j}^{0} + \mathbf{n}_{s} \left ( r, z \right ) \left ( \mathbf{w}_{j}^{\infty } - \mathbf{w}_{j}^{0} \right ), $$

(50)

where $\mathbf{w}_{j}^{\infty }$ and $\mathbf{w}_{j}^{0}$ are the displacement matrices of the $j$th element nodes at $t =\infty$ and $t =0$, respectively, and $\mathbf{n}_{s} \left ( r, z \right )$ represents the shape functions of the elements.

The strain matrix of the $i$th element is defined as follows:

$$ \boldsymbol{\varepsilon }_{i} \left ( r, z, t \right ) = \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{0} + \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} f \left ( t \right ). $$

(51)

Based on Eqs. (50)–(51), the strain rate, strain variations, and displacement variations are defined as follows:

$$ \begin{aligned} \dot{\boldsymbol{\varepsilon }}_{i} \left ( r, z, t \right ) &= \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} \dot{f} \left ( t \right ), \\ \delta\boldsymbol{\varepsilon }_{i} \left ( r, z, t \right ) &= \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} \delta f \left ( t \right ), \\ \delta w_{j} \left (r, z, t \right ) &= \mathbf{n}_{s} \left ( r, z \right ) \mathbf{u}_{j}^{c} \delta f \left ( t \right ). \end{aligned} $$

(52)

We obtain the following equilibrium equation of the circular viscoelastic solid:

$$\begin{aligned} &\sum _{i =1}^{N} 2 \pi \int _{A}^{0} \left ( \mathbf{u}_{i}^{0 T} \mathbf{B}^{T} \mathbf{D} \left ( t \right ) \mathbf{B} \mathbf{u}_{i}^{c} + \int _{0^{+}}^{t} \mathbf{u}_{i}^{c T} \mathbf{B}^{T} \dot{f} \left ( m \right ) \mathbf{D} \left ( t - m \right ) \mathbf{B} \mathbf{u}_{i}^{c} {dm} \right ) r {dA} \\ &\quad = \sum _{j =1}^{N_{s}} 2 \pi \int _{R}^{0} q \mathbf{n}_{s} \mathbf{u}_{j}^{c} r {dr}, \end{aligned}$$

(53)

where $N$ and $N_{s}$ represent the numbers of 2D and 1D elements, respectively.

1.6 A.6 Calculating the stress in the circular viscoelastic prism

As discussed in Sect. 2.5, we write the strain matrix of the $i$th element as follows:

$$ \boldsymbol{\varepsilon }_{i} \left ( r, z, t \right ) = \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{0} + \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} f \left ( t \right ). $$

(54)

By inserting Eq. (54) into Eq. (3) we obtain the equation

$$ \boldsymbol{\sigma }_{i} \left ( r, z, t \right ) = \mathbf{D} \left ( t \right ) \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{0} + \int _{0^{+}}^{t} \mathbf{D} \left ( t - m \right ) \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} \dot{f} \left ( m \right ) {dm}. $$

(55)

Similarly to Sect. 2.7, we obtain the following stress matrix of the $i$th element:

$$ \boldsymbol{\sigma }_{i} \left ( r, z, t \right ) = \mathbf{D} \left ( t \right ) \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{0} + \mathbf{D}_{fc} \mathbf{B} \left ( r, z \right ) \mathbf{u}_{i}^{c} $$

(56)

where

$$ \mathbf{D}_{fc} = \int _{0^{+}}^{t} \mathbf{D} \left ( t - m \right ) \dot{f} \left ( m \right ) dm = K \left [ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \textstyle\begin{array}{c} d_{f 1}\\ d_{f 2}\\ \textstyle\begin{array}{c} 0\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} d_{f 2}\\ d_{f 1}\\ \textstyle\begin{array}{c} 0\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c@{\quad }c} \textstyle\begin{array}{c} 0\\ 0\\ \textstyle\begin{array}{c} d_{f 3}\\ 0 \end{array}\displaystyle \end{array}\displaystyle & \textstyle\begin{array}{c} 0\\ 0\\ \textstyle\begin{array}{c} 0\\ d_{f 3} \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \end{array}\displaystyle \right ]. $$

(57)

Appendix B

The secant method algorithm consists of the following two steps:

Selecting two points in the intended range and calculating the next point by the equation

$$ \alpha _{n +1} = \frac{\alpha _{n -1} \times J \left ( \alpha _{n} \right ) + \alpha _{n} \times J \left ( \alpha _{n -1} \right )}{J \left ( \alpha _{n} \right ) - J \left ( \alpha _{n -1} \right )},\quad n =1,2,3,\dots . $$

(58)

If $J \left ( \alpha _{n} \right ) = J \left ( \alpha _{n -1} \right )$, then $\alpha _{n}$ is considered slightly larger.

$$ \text{If}\ J \left ( \alpha _{n} \right ) = J \left ( \alpha _{n -1} \right ),\quad \text{then}\quad \alpha _{n} = \alpha _{n} + \varepsilon , $$

(59)

where $\varepsilon$ is a small positive number.

Calculations end if $\left \vert J \left ( \alpha _{n} \right ) \right \vert$ is sufficiently close to zero and $\alpha _{n}$ and $\alpha _{n +1}$ are sufficiently close, that is,

$$ \left \vert J \left ( \alpha _{n} \right ) \right \vert < \varepsilon _{1} ,\qquad \left \vert \alpha _{n +1} - \alpha _{n} \right \vert < \varepsilon _{2}, $$

(60)

where $\varepsilon _{1}$ and $\varepsilon _{2}$ are two small positive numbers.

The algorithm returns to Step 2 if the condition in Step 3 is not satisfied.

Appendix C

In this section, we describe the determination of viscoelastic material parameters in the Abacus software (Abaqus 6.14 2014, Documentation-HTML).

In the Abaqus software the time-dependent shear and bulk moduli are formulated by multiplying the elastic moduli at time zero by the dimensionless relaxation functions as follows:

$$\begin{aligned} & G_{R} \left ( t \right ) = G_{0} \left ( 1- \sum _{k=1}^{N} \overline{g}_{k}^{p} \left ( 1- e^{{-t} / {\tau _{k}}} \right ) \right ), \end{aligned}$$

(61)

$$\begin{aligned} &K_{R} \left ( t \right ) = K_{0} \left ( 1- \sum _{k=1}^{N} \overline{k}_{k}^{p} \left ( 1- e^{{-t} / {\tau _{k}}} \right ) \right ). \end{aligned}$$

(62)

Selecting one term from the series results in

$$\begin{aligned} & G_{R} \left ( t \right ) = G_{0} \left ( 1- \overline{g}_{1}^{p} + \overline{g}_{1}^{p} e^{{-t} / {\tau _{1}}} \right ), \end{aligned}$$

(63)

$$\begin{aligned} &K_{R} \left ( t \right ) = K_{0} \left ( 1- \overline{k}_{1}^{p} + \overline{k}_{1}^{p} e^{{-t} / {\tau _{1}}} \right ), \end{aligned}$$

(64)

where $\overline{g}_{1}^{p}$ and $\overline{k}_{1}^{p}$ are the constant parameters of the shear and bulk moduli, respectively, $\tau _{1}$ is the relaxation time, and $G_{0}$ and $K_{0}$ are the shear and bulk moduli at time zero.

So, the nondimensional shear relaxation modulus can be written as

$$ \frac{G_{R} \left ( t \right )}{G_{0}} =1- \overline{g}_{1}^{p} + \overline{g}_{1}^{p} e^{{-t} / {\tau _{1}}}. $$

(65)

In the present study, we assume the bulk modulus to be constant:

$$ K \left ( t \right ) =K. $$

(66)

Also, in the present study, we define the dimensionless shear modulus as follows:

$$ \frac{G \left ( t \right )}{\left ( 3K/2 \right )} =\omega \left ( t \right ) = c_{1} + c_{2} e^{- {t} / {t_{s}}} . $$

(67)

The comparison of Eqs. (65) and (67) yields

$$ \begin{aligned} c_{1} =1- \overline{g}_{1}^{p},\qquad c_{2} &= \overline{g}_{1}^{p},\ t_{s} = \tau _{1},\qquad G_{0} = \frac{3K}{2}, \\ &c_{1} + c_{2} =1 \end{aligned} $$

(68)

and the comparison of Eqs. (64) and (66) yields

$$ \overline{k}_{1}^{p} =0, \qquad K_{0} =K. $$

(69)