Mixed virtual element formulations for incompressible and inextensible problems

Abstract

Locking effects can be a major concern during the numerical modelling of elastic materials, especially for large strains. Those effects arise from volumetric constraints such as incompressibility or anisotropic effects of the underlying material class. One particular solution strategy is to employ mixed formulations, which provide solutions tailored to the specific locking phenomena at hand. While being in general powerful, one drawback of such solution strategies is that the provided strategy to overcome locking is often tied or limited to some specific topological type of finite elements and thus forfeiting generality. Contrary the Virtual Element Method (VEM) benefits by definition from allowing arbitrary element shapes and number of nodes at element level. A variety of approaches for the treatment of locking phenomena for hyperelastic material is formulated in this contribution for the Virtual Element Method with a low order ansatz. Key ingredient to the implementation of multi-field mixed principles in VEM is the consideration of only one constant variable per field and one corresponding Lagrange multiplier over the entire virtual element. Hereby the stabilization contribution utilizes the mixed formulation but shares the element-wise constant variable with the projection part of the virtual element. A direct consequence of this rather simple implementation strategy is the combination of powerful mixed formulations with a computational approach that is able to treat general element shapes. The proposed formulations are tested with regard to structured mesh at standard examples in computational mechanics as well as at specific computational engineering applications where also unstructured meshes are utilized.

A Appendix

1.1 A.1 Computational homogenization: sensitivity analysis for mixed virtual element formulations

The homogenization problem is schematically depicted in Fig. 19. In a multi-scale scheme the constitutive response is modelled at the microscopic length scale while the load is usually applied at the macroscopic length scale. Boundary conditions at microscopic length scales are defined for example via a deformation measure from the macroscopic material point to which the representative volume element (RVE) is attached. Extensive works on such multi-scale scheme formulations in the framework of the $\hbox {FE}^{2}$ are represented in [37, 45, 46, 51] or for non-linear problems in e.g. [35] among many others.

Fig. 19

Schematic illustration of a multi-scale scheme: localization and homogenization of a microscopic computational problem

The unknowns of the non-linear homogenization problem at hands are the macroscopic stress ${\varvec{P}}_{M}$ and its sensitivities ${\mathbb {A}}_{M}=\textrm{D}{\varvec{P}}_{M}/\textrm{D}{\varvec{F}}_{M}$ with repsect to the macroscopic deformation gradient ${\varvec{F}}_{M}$. Those quantities are usually provided by evaluation of the microscopic problem at hands. The homogenization approach is based on the Hill-Mandel condition for large deformations (see i.e. [1]), further assumptions are the absence of body forces at the microscopic domain and perfect bounded grains i.e., no displacement jumps over grain boundaries. The Hill-Mandel condition requires the equivalence of the variation of virtual work, performed at a macroscopic material point ${\varvec{X}}_{M}$ to be equal to the variation of virtual work that is performed in an averaged sense over the entire RVE. In hand with suitable Dirichlet-type boundary conditions, applied to the entire boundary of the RVE $\partial \varvec{\Omega }$ it follows

$$\begin{aligned} \begin{aligned}&{\varvec{P}}_{M}\varvec{:}\delta {\varvec{F}}_{M}\overset{!}{=}\frac{1}{\Vert \varvec{\Omega }\Vert }\int _{\varvec{\Omega }}{\varvec{P}}\varvec{:}\delta {\varvec{F}}\textrm{d}V=\langle {\varvec{P}}\varvec{:}\delta {\varvec{F}}\rangle \\&=\langle {\varvec{P}}\rangle \varvec{:}\delta {\varvec{F}}_{M},\quad \langle (\cdot )\rangle =\frac{1}{V}\int _{\varvec{\Omega }}(\cdot )\textrm{d}V,\\&\Rightarrow {\varvec{P}}_{M}=\langle {\varvec{P}}\rangle . \end{aligned} \end{aligned}$$

(48)

The authors employ linear Dirichlet boundary conditions at the microscale

$$\begin{aligned} {\varvec{u}}({\varvec{X}}_{bc})=({\varvec{F}}_{M} -{\varvec{1}}){\varvec{X}}_{bc},\quad \{{\varvec{X}}_{bc}:{\varvec{X}}\subset {\mathbb {R}}^{d}\wedge {\varvec{X}}\in \partial \varvec{\Omega }\}, \end{aligned}$$

(49)

with d being the spatial dimension of the problem at hands. Combination of (48) and (49) yields

$$\begin{aligned} \begin{aligned}&{\varvec{P}}={\varvec{P}}({\varvec{F}}({\varvec{u}}({\varvec{F}}_{M}))),\\&\quad {\mathbb {A}}_{M}=\frac{\textrm{D}{\varvec{P}}_{M}}{\textrm{D}{\varvec{F}}_{M}}=\frac{\textrm{D}\langle {\varvec{P}}\rangle }{\textrm{D}{\varvec{F}}_{M}}=\frac{1}{\Vert \varvec{\Omega }\Vert }\displaystyle \int _{\varvec{\Omega }}\frac{\textrm{D}{\varvec{P}}}{\textrm{D}{\varvec{F}}_{M}}\textrm{d}V,\\&\Rightarrow \frac{\textrm{D}{\varvec{P}}}{\textrm{D}{\varvec{F}}_{M}}=\frac{\partial {\varvec{P}}}{\partial {\varvec{F}}}\frac{\partial {\varvec{F}}}{\partial {\varvec{u}}}\frac{\textrm{D}{\varvec{u}}}{\textrm{D}{\varvec{F}}_{M}}. \end{aligned} \end{aligned}$$

(50)

From (50) it is obvious that the only unknown term for the derivation of ${\mathbb {A}}_{M}$ is the sensitivity $\textrm{D}{\varvec{u}}/\textrm{D}{\varvec{F}}_{M}$ of microscopic displacements with respect to the macroscopic deformation gradient. The sensitivity analysis is based on formulations from [29, 55, 66] and in the following subsequently extended to the presented virtual element scheme with a mixed formulation. States of displacement ${\varvec{u}}$ drive the mechanical problem at hands and being the nodal unknowns ${\textbf{p}}$ in the virtual element formulation. Starting from (38) (resp. (39)), the solution to the primary problem is found when ${\textbf{R}}={\varvec{0}}$. Rewriting the residual leads to

(51)

Hereby $\varvec{\Phi }$ denotes the set of design (sensitivity) parameters, being the individual components of ${\varvec{F}}_{M}$ in the constitutive framework under consideration. Note, that ${\textbf{R}}_{int}$ exhibits a dependency on the condensation variables, stored in ${\textbf{h}}_{s}$. These variables are considered as history-variables in the following development, furthermore they have a dependency on the design parameters $\varvec{\Phi }$ which results in additional terms for the sensitivity analysis. Application of the total derivative to (51) leads to

(52)

Note, that the bracket-term in (52) yields exactly the factorized tangent matrix before condensation from (38). ${\mathbb {R}}_{t}$ denotes a first order pseudo load term, $\textrm{D}{\varvec{u}}_{bc}/\textrm{D}\varvec{\Phi }$ is the sensitivity of constrained nodal displacements with respect to the design paramters and has to be given as an input. Consideration of boundary conditions from (49) leads to a rather simple solution $\textrm{D}{\varvec{u}}_{bc}/\textrm{D}\varvec{\Phi }={\varvec{X}}_{bc}$. (52) yields a system of algebraic equations for the unknown sensitivities $\textrm{D}{\varvec{u}}_{h}/\textrm{D}{\varvec{F}}_{M}$ that is solved in a straight forward manner by a Newton–Raphson procedure. From here, ${\mathbb {A}}_{M}$ can be computed. It is worth mentioning, that the presented scheme restores an algorithmic consistent linearized expression for ${\mathbb {A}}_{M}$. Furthermore the authors want to note that the presented sensitivity analysis is a rather computational efficient scheme since the solution of the sensitivity problem is sought after the primal problem was solved, which provides the already inverted and factorized tangent matrix ${\textbf{K}}$ from the last (and successful) solution step and thus only ${\mathbb {R}}_{t}$ has to be built. For more details, the interested reader is forwarded to [31] and references therein. The effective algorithmic consistent linearized material matrix for VEM and FEM approaches is computed by

$$\begin{aligned} \begin{aligned}&\mathrm{VEM:}{\mathbb {A}}_{h}=\frac{\textrm{D}\langle {\varvec{P}}\rangle }{\textrm{D}{\varvec{F}}_{M}}\\&\quad =\frac{1}{\Vert \varvec{\Omega }\Vert }\sum _{E\in \varvec{\Omega }} \left( (1-\beta |_{E})V_{E}\frac{\partial {\varvec{P}}_{\pi }|_{E}}{\partial {\varvec{F}}_{\pi }|_{E}} \frac{\partial {\varvec{F}}_{\pi }|_{E}}{\partial {\varvec{u}}_{{h}|_{E}}}\right. \\&\quad \left. +\beta |_{E}\sum _{{\mathfrak {T}}\in E}w_{{\mathfrak {T}}}J_{{\mathfrak {T}}}\frac{\partial {\varvec{P}}_{{\mathfrak {T}}} |_{E}}{\partial {\varvec{F}}_{{\mathfrak {T}}}|_{E}}\frac{\partial {\varvec{F}} _{{\mathfrak {T}}}|_{E}}{\partial {\varvec{u}}_{{h}|_{E}}}\right) \frac{\textrm{D}{\varvec{u}}_{{h}|_{E}}}{\textrm{D}{\varvec{F}}_{M}},\\&\mathrm{FEM:}{\mathbb {A}}_{h}=\frac{\textrm{D}\langle {\varvec{P}}\rangle }{\textrm{D}{\varvec{F}}_{M}}\\&\quad =\frac{1}{\Vert \varvec{\Omega }\Vert }\sum _{E\in \varvec{\Omega }}\left( \sum _{g_{p}\in E}w_{g_{p}}J_{g_{p}}\frac{\partial {\varvec{P}}|_{g_{p}}}{\partial {\varvec{F}}|_{g_{p}}} \frac{\partial {\varvec{F}}|_{g_{p}}}{\partial {\varvec{u}}_{{h}|_{E}}}\right) \frac{\textrm{D}{\varvec{u}}_{{h}|_{E}}}{\textrm{D}{\varvec{F}}_{M}}. \end{aligned} \end{aligned}$$

(53)

$w_{{\mathfrak {T}}}=w_{g_{p}}$ denotes the integration point weight of the integration point $g_{p}\in {\mathfrak {T}}$, $J_{{\mathfrak {T}}}=J_{g_{p}}$ is the jacobian at integration point $g_{p}\in {\mathfrak {T}}$. Note, that linear triangles/tetrahedra are used for the stabilization mesh, thus $n_{g_{p}}=1$.

1.2 A.2 Additional results

For the 2D Cook’s membrane problem see Fig. 20. Additional results for the coloured plots of $[{\mathbb {A}}_{h}|_{g}]_{ij}$ are depicted in Figs. 21, 22 and 23.

Fig. 20

Results for: a, c, e the convergence and b, d, f the error of the 2D Cook’s membrane problem subjected to the Neo-Hookean transversely isotropic material response with different degrees of anisotropy in (11) (resp. (18))

Fig. 21

Grain-averaged predictions $[{\mathbb {A}}_{h}|_{g}]_{2j}=\textrm{D}[\langle {\varvec{P}}\rangle _{kl}|_{g}]_{j}/\textrm{D}\varvec{\Phi }_{2}$ of the in-plane entries $j=(2,3,4)$: VEM VO-CoTr-JP (a, d, g), benchmark FEM T2 fine (b, e, h), grain-wise absolute deviation $\Delta [{\mathbb {A}}_{h}|_{g}]_{2j}$ (c,f,i)

Fig. 22

Grain-averaged predictions $[{\mathbb {A}}_{h}|_{g}]_{ij}=\textrm{D}[\langle {\varvec{P}}\rangle _{kl}|_{g}]_{j}/\textrm{D}\varvec{\Phi }_{i}$ of the in-plane entries $i=(3,4),\;j=(3,4)$: VEM VO-CoTr-JP (a, d, g), benchmark FEM T2 fine (b,e,h), grain-wise absolute deviation $\Delta [{\mathbb {A}}_{h}|_{g}]_{ij}$ (c, f, i)

Fig. 23

Grain-averaged predictions $[{\mathbb {A}}_{h}|_{g}]_{i5}=\textrm{D}[\langle {\varvec{P}}\rangle _{kl}|_{g}]_{5}/\textrm{D}\varvec{\Phi }_{i}$ of the out-of-plane entries: VEM VO-CoTr-JP (a, d, g, j), benchmark FEM T2 fine (b, e, h, k), grain-wise absolute deviation $\Delta [{\mathbb {A}}_{h}|_{g}]_{i5}$ (c, f, i, l)