documents

Overview

Many rubber-like materials are not only near-incompressible in nature, but also exhibit a viscoelastic response (within the tested load and time scales). In this example, we extend the near-incompressible rate-independent constitutive used in step-44 (which implements three-field quasi-static quasi-incompressible finite elasticity) to one that displays rate-dependent behavior. It may seem that there is a contradiction of terms here, so lets clarify by saying the the problem remains "quasi-static" in the sense that inertial terms remain insignificant, even though the material response itself is rate-dependent. This implies that, for these fictitious material parameters, it is assumed that the timescale for material relaxation is much longer than that of elastic wave propagation.

We've also taken the opportunity to extend the code first shown in step-44 to parallel (the primary reason for this contribution), using Metis as a grid partitioner and Trilinos for linear algebra. As a motivation as to why one might still choose to use Metis (also associated with parallel::shared::Triangulation in step-18, although this triangulation is not used in this instance) over p4est (also associated with parallel::distributed::Triangulation) as a grid partitioner, at this point in time it is not possible to use the hp finite-element in conjunction with the distributed grid, meaning that this code could not, for example, be readily extended to the application shown in

The discerning reader will observe that we've chosen to employ deal.II's built in solvers as opposed to using Trilinos solvers. This is because the system matrices K_Jp and K_pJ, although block diagonal and well conditioned, and for some reason (perhaps pertaining to the negative definite nature of these blocks, or that the entries are very small in magnitude) Trilinos solvers are not sufficiently robust to compute inverse matrix-vector multiplication with. We do stress, however, that to date no great attempt has been made by the author to overcome this issue other than by making an entirely different choice of solver.

Finite deformation of a thin strip with a hole.

Various permutations of the problem of an elastomeric strip with a centered cut-out can be found in the literature for solid mechanics, in particular (but not limited to) that pertaining to

Here we implement another permutation (one that is not necessarily benchmarked elsewhere), simply for demonstration purposes. The basic problem configuration is summarized in the following image.

Problem geometry

A thin strip of material with a circular hole is (in 3d) constrained in the Z direction and loaded in the direction of its long edge. In our implementation, this load is applied to the +Y surface and may either be displacement control (a Dirichlet condition) or a traction load (a Neumann boundary condition). Due to the symmetry of both the geometry and load, the problem can be simplified by modeling only an eighth of the geometry in 3d or a quarter in 2d. By doing so, it it necessary to then implement symmetry conditions on the surfaces coincident with the X-Z and Y-Z planes (and the X-Y plane in 3d). The +X surface, and that of the hole itself, remain traction-free.

In three dimensions, the geometry (and a potential partitioning over 8 processors) looks as follows:

3d grid with partitioning

Note that, for this particular formulation, the two-dimensional case corresponds to neither plane-strain nor plane-stress conditions.

Requirements

Compiling and running

Similar to the example programs, run

cmake -DDEAL_II_DIR=/path/to/deal.II .

in this directory to configure the problem.
You can switch between debug and release mode by calling either

make debug

or

make release

The problem may then be run in serial mode with

make run

and in parallel (in this case, on 4 processors) with

mpirun -np 4 ./viscoelastic_strip_with_hole

This program can be run in 2d or 3d; this choice can be made by making the appropriate changes in the main() function.

Reference for this work

If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation:

Acknowledgements

The support of this work by the European Research Council (ERC) through the Advanced Grant 289049 MOCOPOLY is gratefully acknowledged by the author.

The derivation of the finite-element problem, namely the definition and linearization of the residual and their subsequent discretization are quite lengthy and involved. However, this is already detailed in step-44, some of the aforementioned literature, as well as

and need not be repeated here. As for the viscoelastic constitutive law (which satisfies the dissipation inequality through the definition of an internal variable that we denote as Q in the code), this is derived and presented in detail in

In particular, the relevant equations from Linder et al.'s work that are implemented in this work are equations 47, 54 and 56. Note that the time discretization for the rate-dependent internal variable for this single dissipative mechanism is only first-order accurate.

Results

To begin, here is a comparison of the initial grid for the 2d version of the problem (mirrored about the x-axis and rotated 90 degrees anti-clockwise)

Initial grid

and of the final, displaced grid after the load has been applied and the material is in a (near-completely) relaxed state.

Displaced grid

These results, as well as those that follow, were produced using the following material properties:

while the boundary conditions were configured in the following manner:

The following chart highlights the displacement of several vertices and clearly illustrates the viscoelastic nature of the material.

Problem geometry

During the initial phase, the load is applied over a period of time much shorter than the material's characteristic relaxation time. The material therefore exhibits a very stiff response, and relaxes as the load remains constant for the remainder of the simulation. This deformation that occurs under constant load is commonly known as creep.

We've been lazy and stopped the simulation slightly prematurely, but it is clear that the material's displacement is moving asymptotically towards a equilibrium solution. You can also check what the true resting load state is by removing the dissipative mechanism (setting its shear modulus to zero) and rerunning the simulation with the material then exhibiting rate-independent behavior. Note that in this case, the number of time step over which the load is applied must likely be increased in order to ensure stability of the nonlinear problem.

Animations

Just in case you found the presented chart a little dry and undigestible, below are a couple of animations that demonstrate the viscoelastic nature of the material in a more visually appealing manner.

Animation of the evolution of the displacement field Displacement field

Animation of the evolution of the pressure field Pressure field