A Unified Determinant-Preserving Formulation for Compressible/Incompressible Finite Viscoelasticity

Abstract

This paper presents a formulation alongside a numerical solution algorithm to describe the mechanical response of bodies made of a large class of viscoelastic materials undergoing arbitrary quasistatic finite deformations. With the objective of having a unified formulation that applies to a wide range of highly compressible, nearly incompressible, and fully incompressible soft organic materials in a numerically tractable manner, the viscoelasticity is described within a Lagrangian setting by a two-potential mixed formulation. In this formulation, the deformation field, a pressure field that ensues from a Legendre transform, and an internal variable of state $F^v$ that describes the viscous part of the deformation are the independent fields. Consistent with the experimental evidence that viscous deformation is a volume-preserving process, the internal variable $F^v$ is required to satisfy the constraint det $F^v=1$. To solve the resulting initial-boundary-value problem, a numerical solution algorithm is proposed that is based on a finite-element (FE) discretization of space and a finite-difference discretization of time. Specifically, a Variational Multiscale FE method is employed that allows for an arbitrary combination of shape functions for the deformation and pressure fields. To deal with the challenging non-convex constraint det $F^v=1$, a new time integration scheme is introduced that allows to convert any explicit or implicit scheme of choice into a stable scheme that preserves the constraint det $F^v=1$ identically. A series of test cases is presented that showcase the capabilities of the proposed formulation.

Keywords: Elastomers; Finite deformations; Stabilized finite elements; Stable ODE solvers.

Figure 1:

Rheological representation of the two-potential model (6)-(7) for viscoelasticity.

Figure 2:

Case 1: Error convergence of (a) bE and (b) the proposed modified bE.

Figure 3:

Case 1: Error convergence of (a) RK5 and (b) the proposed modified RK5.

Figure 4:

Case 2: Error convergence of (a) bE and (b) the proposed modified bE.

Figure 5:

Case 2: Error convergence of (a) RK5 and (b) the proposed modified RK5.

Figure 6:

Evolution in time t of $\text{det}\hspace{0.17em}{\mathbf{Y}}^b$ and $\text{det}\hspace{0.17em}\mathbf{Y}$ for (a) Case 1 and (b) Case 2.

Figure 7:

Error convergence in the measure $\left|\text{det}\hspace{0.17em}{\mathbf{Y}}^b-1\right|$ of (a) fE and (b) bE. Evolution in time $t$ of $\text{det}\hspace{0.17em}{\mathbf{Y}}^b$ and $\text{det}\hspace{0.17em}\mathbf{Y}$ for (c) fE and (d) bE.

Figure 8:

Comparisons with the results in [25] for the response of VHB 4910 under: (a) uniaxial tension loading/unloading at various constant stretch rates and (b) a single-step relaxation with two different applied stretches.

Figure 9:

Schematic of the indentation of a block made of viscoelastic elastomers with disparate compressibilities.

Figure 10:

Contour plots of the displacement field $u_3\left(\mathbf{X},t\right)$ at $t=T=2$ s for T4 elements and several meshes of increasing refinement.

Figure 11:

Contour plots of the displacement field $u_3\left(\mathbf{X},t\right)$ at $t=T=2$ s for H8 elements and several meshes of increasing refinement.

Figure 12:

Evolution in time $t$ of $\text{det}\hspace{0.17em}\mathbf{C}$ averaged over an octant of the cube directly below the applied load for T4 elements and the proposed time integration scheme with base scheme (a) bE and (b) RK5.

Figure 13:

Evolution in time $t$ of $\text{det}\hspace{0.17em}\mathbf{C}$ averaged over an octant of the cube directly below the applied load for H8 elements and the proposed time integration scheme with base scheme (a) bE and (b) RK5.

Figure 14:

Contour plots of $\text{det}\hspace{0.17em}{\mathbf{C}}^v$ over the deformed configuration at the final instance $t=T=2$ s of the applied loading for T4 elements. The results are shown for (a) the bE scheme and (b) the proposed modified bE scheme.

Figure 15:

Contour plots of $\text{det}\hspace{0.17em}{\mathbf{C}}^v$ over the deformed configuration at the final instance $t=T=2$ s of the applied loading for H8 elements. The results are shown for (a) the bE scheme and (b) the proposed modified bE scheme.

Figure 16:

Evolution in time of $\text{det}\hspace{0.17em}{\mathbf{C}}^v$ at the integration point closet to $\left(X_1,X_2,X_3\right)=\left(0,0,0.5\right)\text{mm}$. The results are shown for (a) the bE and the proposed modified bE schemes and for (b) the RK5 and the proposed modified RK5 schemes.

Figure 17:

Geometry of the specimen, mesh, and contour plot of the first Piola-Kirchoff stress $P_{33}\left(\mathbf{X},t\right)$ over the deformed configuration at the final time $t=T=50$ s of the applied loading.

Figure 18:

Comparison of the force-time response predicted by the simulation with the indentation experimental data and simulation reported in [1] for porcine liver.

Figure 19:

Specimen containing 200 randomly distributed initially spherical bubbles of monodisperse size at volume fraction $c=0.15$ and its FE discretization with approximately 1.4 million linear tetrahedral (T4) elements.

Figure 20:

Response of the suspension of bubbles under simple shear. The plots show the resulting average Cauchy shear stress ${\overline{T}}_{12}$ as a function of the applied Eulerian shear rate ${\overline{D}}_{12}$. For direct comparison, the plot includes the corresponding response (dashed line) of the underlying elastomer without the bubbles.