Large strain Flory’s decomposition for Lagrangian modeling of viscoleastic solids and compressive fluids

Abstract

The fundamental difference in the solution of solids and fluids relies on the respective constitutive laws. Based on the Rivlin-Saunders-Düster-Hartmann hyperelastic model and using the Flory’s strain decomposition, we present a new total Lagrangian viscoelastic constitutive model for both Kelvin-Voigt viscoelastic solids and free-surface compressive viscous isothermal fluids. A dissipative viscous virtual work is written as a function of the time rate of isochoric invariants and its relation with the viscous stress is derived. Local time derivatives are solved by backward finite difference, allowing a consistent tangent viscoelastic constitutive tensor. The virtual work principle is used to write the weak equilibrium equation and its position-based finite element counterpart. Dynamic time integration is carried out by the Newmark $\beta$ method and the Newton-Raphson procedure is used to solve time steps. The formulation is validated against experimental and numerical literature results revealing good precision. Additional examples are shown in order to demonstrate the applicability and future possibilities of the technique.

Keywords
Unified solid-fluid model; Kelvin-Voigt visco-hyperelasticity; Flory’s decomposition; Positional FEM; Total Lagrangian formulation

Graphical Abstract

1 INTRODUCTION

Regarding solid applications, the development of the finite element method (FEM) is quite old and, among others, pioneers works should be mentioned and briefly commented: Courant (1942) introduced the idea of the minimization of a functional using linear approximation over sub-regions, Argyris (1954) developed the matrix theory of structures for the discrete elements leading to concepts as stiffness and flexibility, Turner et al. (1956) started discussing global stiffness of truss elements and introduced triangular plate elements, Clough (1960) introduced the name Finite Element Method in an important work on dynamics, Clough (1989), and Zienkiewicz and Cheung (1965) extended the FEM for no-structural applications and torsion. Originally, the solution of linear structural problems, with small displacements and simple constitutive relations, has been the focus of FEM researches. Formulations considering large displacements and small strains are part of FEM evolution regarding solids and structures, for a brief description one may consult Bischoff and Ramm (2000), Sansour and Bednarczyk (1995), Gruttmann and Wagner (2001) and Benson et al. (2011). These works were interested in solving thin plates and shells developing large displacements and rotations.

More challenging problems, involving nonlinear constitutive laws as plasticity and damage mechanics, have been widely studied using FEM, even though considering small displacements and strains (Havner (1966), Miehe et al. (2013)). However, considering small strains avoids modeling hyperelastic materials as polymers and problems related to manufacturing and general structural collapse as, for example, cold forming and crashworthiness.

Problems involving large strains, large displacements and hyperelastic nonlinear constitutive relations, also solved by FEM, bring important contributions to the understanding of highly deformable solids and general applications as one can see, for example, in Holzapfel and Simo (1996) for isotropic finite strain elasticity, Gasser and Holzapfel GA (2002) for orthotropic and biological elastic applications, Vergori et al. (2013) and Latorre and Montans (2015) for anisotropic finite elasticity. As the constitutive equations of hyperelastic materials will be the focus of our derivations, it is important to mention that its constitutive models can be classified mainly in 3 groups: (i) phenomenological, which are based on experimental observations as the Mooney-Rivlin model (Mooney (1940), Rivlin (1948a) and Rivlin (1948b)) and the models proposed by Yeoh (1990) and Ogden (1972); (ii) mechanistic, which considers the structure of the mater as the Arruda and Boyce (1993) model and the Neo-Hookean model (Treloar (1943)); and (iii) hybrids, that mix both groups, see Gent (2001) for details. One can mention viscoelastic model, where strain rates are also considered (Argyris et al. (1991), Chen et al. (1993), Holzapfel (1996), Shutov et al. (2013) and Pascon and Coda (2017)). When large strains and nonlinear constitutive laws (plasticity or viscoplasticity) are considered, the idea of near incompressibility and flux laws appear.

In general, formulations that consider finite strain plasticity use local plastic potentials to define flux direction and the Kröner-Lee (Kröner (1959) and Lee (1969)) multiplicative decomposition to split elastic and plastic strain. This strategy results in difficult expressions and in an intermediary plastic space (Simo (1992), Reese and Wriggers (1997), Jiao and Fish (2018) and Hudobivnik et al. (2019)). Some interesting FEM applications using plastic potentials, as metal forming and stamping, can be seen for example in Chung et al. (1998) and Schwarze et al. (2011). The same Kröner-Lee multiplicative decomposition is also used to describe finite viscoelasticity or viscoplasticity (Pascon and Coda (2015) and Pascon and Coda (2017)). Although considering large strain viscoplasticity, none of the above mentioned works extended their formulation to model fluids, as it will be done originally here. At the same time the evolution of FEM modeling of fluids, briefly described in what follows, does not goes in the direction of a unification with solids modeling.

Regarding fluid applications, a general overview of several computational frameworks can be found in textbooks as Anderson (1995), Chung (2002), Zienkiewicz et al. (2005) and Reddy and Gartling (2010). Due to some aspects, such as the facility to use arbitrarily unstructured meshes and, particularly, due to the simplicity of boundary conditions enforcement over complex boundaries, the FEM has been conquering more space in fluid mechanics. However, fluid dynamics analyses under Eulerian description may present dominant convection, so the direct application of standard Galerkin method, as the FEM, to solve the governing equations leads to nonsymmetrical matrices with spurious behavior in solutions (Brooks and Hughes (1982) and Elias et al. (2006)). Many researchers have developed techniques to increase FEM stability and robustness for fluid dynamics, such as Tezduyar (1992), Tezduyar and Senga (2006), Akkerman et al. (2008), among others.

For free surface flow problems it is not possible to apply a pure Eulerian description, as the fluid domain is dynamically deformed, denoting a moving-boundary problem. In this context, Arbitrary Lagrangian Eulerian (ALE) stabilized formulations have been developed by Donea et al. (1982), Franci et al. (2016) and Duarte et al. (2004); and the stabilized space-time formulation by Tezduyar et al. (1992). Such methods allow deforming the fluid domain discretization by mesh updating. However, this mesh updating respects the fluid movement only at the boundary and fluid flow is allowed between general internal elements, eliminating any possibility of solid modeling.

As an alternative to ALE, the updated Lagrangian description has been applied together with particle methods concepts to develop the Particle Finite Element Method (PFEM) (Idelsohn et al. (2004), Idelsohn et al. (2006) and Idelsohn et al. (2008)), providing a good solution for free surface flows with topological changes in the fluid domain. This strategy also demands remeshing and needs special attention in defining physical properties. The PFEM is dedicated to model fluids, the remeshing is specialized for linear triangle elements and - during this process - elements are eliminated and/or created, i.e., the mass balance is not preserved.

In this study we propose an original link between solid mechanics and free surface flow fluid mechanics. We take advantage of the finite strain solid mechanics for hyperelastic materials and its tensor algebra, Flory’s decomposition (Flory (1961)), to propose a total Lagrangian formulation covering both Kelvin-Voigt-like viscoelastic solids and viscous fluids. When simulating fluids, the proposed formulation is limited to isothermal-compressible-viscous free-surface flows with finite distortions and free of topological changes (no surface break and no fluid-fluid surface contact). The resulting constitutive law is implemented in a total Lagrangian positional FEM code integrally developed in the Structural Engineering Department of the University of São Paulo. As it is a solid-based computational code, it is important to mention that the local viscous time derivatives are solved by backward finite difference, making possible to write a consistent tangent viscoelastic constitutive tensor. Time integration is carried out by the Newmark $\beta$ method and the Newton-Raphson procedure is used to implicitly solve time steps. It is important to mention that the positional FEM has been developed and published originally in Coda (2003) and applied in various structural subjects as: flexible multi-body dynamics in Coda and Greco (2006), geometrical non linear analysis of shells in Coda and Paccola (2008), inflatable structures in Coda (2009), 3D frames in Coda and Paccola (2011), inverse problems in Coda (2015), in aero space structures in Siqueira and Coda (2019) and finite elastoplasticity in Coda (2021). There are three main advantages of using positional FEM, in order of importance: (i) the use of unconstrained vectors to substitute finite rotations in 3D frame and shell analyses, (ii) the natural presence of a numerical chain rule that facilitates the mapping of geometrical nonlinearities and (iii) as the resulting formulation is total Lagrangian the mass matrix is constant and the time integration can be done in the same way of geometrical linear dynamic analysis.

This work is organized as follows: Item 2.1 describes the weak form of the equilibrium equations, valid for both solid and fluids, detaching the total stress that will be considered for both solid and fluid constitutive models. In item 2.2 we describe the adopted hyperelastic constitutive equations emphasizing the uncoupled strain directions. In item 2.3 the dissipative viscous virtual work is written regarding the time rate of isochoric invariants and both the viscous stress and the tangent viscoelastic constitutive tensor are derived. Section 3 describes the steps to implement the positional FEM computational code. In section 4 low-compressible flows and Kelvin-Voigt-like visco-hyperelastic examples are used to validate the proposed model and to illustrate its possibilities. Acknowledgements and conclusions are given in sections 5 and 6, respectively.

2 DEVELOPMENTS – CONSTITUTIVE MODEL CONCEPTUALIZATION

This section introduces the original Lagrangian finite strain viscous model to be employed in solids and fluids applications. In order to do so, we divide the developments into 3 topics, the first introduces the Lagrangian weak equilibrium (movement) equations detaching its continuum mechanics characteristics, the second introduces the meaning of uncoupled (volumetric and isochorics) Lagrangian strain directions and the third topic presents the viscous constitutive model we seek.

2.1 Equilibrium equations – Strong and weak forms

Starting from the Eulerian local (strong) equilibrium equations, we write the weak form of the equilibrium equations in its Lagrangian version and in a suitable shape for modeling continuum problems by the positional FEM.

The movement equations in the Eulerian description, see for instance Ogden (1984), are written using the Cauchy stress tensor as:

$$\nabla \cdot {\sigma }^T+\overrightarrow{b}=\rho \overrightarrow{\ddot{y}}$$ (1)

and

$$\sigma ={\sigma }^T$$ (2)

in which $y_i$ is a position of a particle inside the continuum, ${\sigma }_{ij}$ is the Cauchy stress component acting on plane $i$ in the direction $j$, $\rho$ is the mass density, ${\ddot{y}}_i$ is the particle acceleration following direction $i$ (material derivative of velocity) and $b_i$ is the $i-th$ body force component. Eq. (2) corresponds to the 3 equilibrium equations in moment that are satisfied for static and dynamic problems by the symmetry of the stress tensor (Ogden(1984)). Eq. (1) contains the 3 translation equilibrium equations in its strong (local) form.

There are many ways to achieve the weak form of Eq. (1), preparing it to FEM implementation. Here, we apply the virtual work concept by the following scalar product,

$$\delta w=\left(\nabla \cdot {\sigma }^T+\overrightarrow{b}-\rho \overrightarrow{\ddot{y}}\right)\cdot \delta \overrightarrow{y}=0$$ (3)

in which $\delta y_i$ is a variation of positions $y_i$.

Integrating Eq. (3) over the spatial domain, the virtual work is written as:

$$\delta W={\displaystyle {\int }_V^{}\left(\nabla \cdot {\sigma }^T+\overrightarrow{b}-\rho \overrightarrow{\ddot{y}}\right)\cdot \delta \overrightarrow{y}dV}=0$$ (4)

and, by distributing the integral, one gets:

$${\displaystyle {\int }_V^{}\rho \overrightarrow{\ddot{y}}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_V^{}\overrightarrow{b}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_V^{}\left(\nabla \cdot {\sigma }^T\right)\cdot \delta \overrightarrow{y}dV=0}$$ (5)

In order to deal with the last term of Eq. (5), we use the following property:

$$\left(\nabla \cdot {\sigma }^T\right)\cdot \delta \overrightarrow{y}=\nabla \left({\sigma }^T\cdot \delta \overrightarrow{y}\right)-{\sigma }^T:\nabla \left(\delta \overrightarrow{y}\right)$$ (6)

Replacing Eq. (6) in the last term of Eq. (5), and applying the divergence theorem, one writes:

$$\delta W={\displaystyle {\int }_V^{}\rho \overrightarrow{\ddot{y}}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_V^{}\overrightarrow{b}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_S^{}{\sigma }^T\cdot \delta \overrightarrow{y}\cdot \overrightarrow{n}dS+{\displaystyle {\int }_V^{}{\sigma }^T:\nabla (\delta \overrightarrow{y})dV}=0}$$ (7)

where $n_j$ is the $j-th$ component of the normal unity vector to the current surface $S$ of the analyzed continuum. Recalling the Cauchy formula, $h_i={\sigma }_{ji}{n}_j$ (in which $h_i$ is the distributed force over the Neumann boundary), and taking into account the stress tensor symmetry, we write:

$$\delta W={\displaystyle {\int }_V^{}\rho \overrightarrow{\ddot{y}}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_V^{}\overrightarrow{b}\cdot \delta \overrightarrow{y}dV}-{\displaystyle {\int }_S^{}\overrightarrow{h}\cdot \delta \overrightarrow{y}dS+{\displaystyle {\int }_V^{}\sigma :\delta \epsilon dV}=0}$$ (8)

with $\delta {\epsilon }_{ij}=(\delta y_{i,j}+\delta y_{j,i})/2$ being a variation of real strain (measured in Eulerian reference), which is mostly known in its rate version ${\dot{\epsilon }}_{ij}=({\dot{y}}_{i,j}+{\dot{y}}_{j,i})/2$ (Ogden(1984) ).Grouping terms of Eq. (8) one has:

$$\delta W=\delta K+\delta P+\delta \Psi$$ (9)

in which $K$ is the kinetic energy, $P$ is the external force potential and $\Psi$ is the Helmholtz free energy (Ogden(1984)) for isothermal states that includes, in this study, the specific strain energy and viscous dissipation. Moreover, Eq. (8) is the weak form of the equilibrium equation for the continuum media in the Eulerian reference.

To obtain the Lagrangian description of Equation (8), the continuity theorem is applied to the first term of Eq. (8),

$${\displaystyle {\int }_V^{}\rho \overrightarrow{\ddot{y}}.\delta \overrightarrow{y}dV}={\displaystyle {\int }_{V_0}^{}{\rho }_0\overrightarrow{\ddot{y}}\cdot \delta \overrightarrow{y}d{V}_0}$$ (10)

and volume forces are considered proportional to density,

$${\displaystyle {\int }_V^{}\overrightarrow{b}\cdot \delta \overrightarrow{y}dV}={\displaystyle {\int }_V^{}\rho \overrightarrow{g}\cdot \delta \overrightarrow{y}dV}={\displaystyle {\int }_{V_0}^{}{\rho }_0\overrightarrow{g}\cdot \delta \overrightarrow{y}d{V}_0}$$ (11)

In equations (10) and (11) we use the well known property $\rho ={\rho }_0/J$ in which $J=\text{d}V/\text{d}{V}_0$ is the determinant of the deformation gradient defined in what follows. In this work, we consider that tractions over Neumann boundary are conservative forces. Cauchy stress tensor $\sigma$, second Piola-Kirchhoff stress tensor $S$, real strain variation $\delta \epsilon$ and Green strain variation $\delta E$ (energetically conjugate to the second Piola-Kirchhoff stress), are related by (Ogden(1984) and Bonet and Wood (1997)):

$$\sigma =\frac{1}{J}AS{A}^T$$ and $$\sigma :\delta \epsilon =S:\delta E$$ (12)

in which $A_{ij}$ is the deformation gradient as usually given in finite strain elasticity, $J$ is its determinant and the Green strain is defined as $E=(A^{t}A-I)/2=(C-I)/2$ with $C=A^{T}A$ being the Cauchy-Green stretch.

From Equations (11) and (12), the week form of the equilibrium equation in the Lagrangian description is written as:

$$\delta W={\displaystyle {\int }_{V_0}^{}{\rho }_0\overrightarrow{\ddot{y}}\cdot \delta \overrightarrow{y}d{V}_0}-{\displaystyle {\int }_{V_0}^{}{\overrightarrow{b}}_{}^0\cdot \delta \overrightarrow{y}d{V}_0}-{\displaystyle {\int }_{S_0}^{}\overrightarrow{h}\cdot \delta \overrightarrow{y}d{S}_0}+{\displaystyle {\int }_{V_0}^{}S:\delta Ed{V}_0=0}$$ (13)

Eq. (13) will be recalled in section 3 to define the FEM formulation, in this study using positions as parameters. It is interesting to advance that the proposed viscous modeling will compose the total stress $S$ together with the hyperelastic stress.

2.2 Hyperelastic constitutive law and Lagrangian uncoupled strain directions

In this topic we describe the multiplicative decomposition of the deformation gradient due to Flory (1961) and, from its result; we introduce the Lagrangian uncoupled strain directions that allows the definition of the proposed viscous model. We also state the adopted hyperelastic part of the complete model.

Cauchy-Green stretch decomposition

Adopting the Green strain as reference, one writes the variation of the Helmholtz free energy of Eq. (9) as:

$$\delta \psi =\frac{\partial \psi }{\partial E}:\delta E$$ with $$\frac{\partial \psi }{\partial E}=S$$ (14)

In hyperasticity it is a current procedure to split the Helmholtz free energy (strain energy in this topic) in two parts, volumetric and isochoric, as follows:

$$\psi ={\psi }^{vol}(\widehat{C})+{\psi }^{iso1}\left({\overline{I}}^1(\overline{C})\right)+{\psi }^{iso2}\left({\overline{I}}^2(\overline{C})\right)$$ (15)

whose parameters are the volumetric $\widehat{C}$ and isochoric $\overline{C}$ parts of the Cauchy-Green stretch tensor $C$. It is worth noting that ${\psi }^{iso1}$ and ${\psi }^{iso2}$ are written, respectively, as a function of the first and second invariants of the isochoric part of the Cauchy-Green stretch tensor. The volumetric $\widehat{C}$ and isochoric $\overline{C}$ parts of the Cauchy-Green stretch tensor are achieved applying the Flory’s decomposition over the deformation gradient $A$ (see Eq. (12)), as follows:

$$A=\widehat{A}\overline{A}$$ (16)

with

$$\widehat{A}=J^{1/3}I\Rightarrow \det \left(\widehat{A}\right)=J$$ (17)

$$\overline{A}=J^{-1/3}A\Rightarrow \det \left(\overline{A}\right)=1$$ (18)

Using the shape of Eq. (16) and taking advantage of the definitions given by Eqs. (17) and (18), we calculate the Cauchy-Green stretch as:

$$C=A^{T}A=J^{2/3}{\overline{A}}^T\overline{A}=J^{2/3}\overline{C}$$ or inversely $$\overline{A}=J^{-2/3}C$$ (19)

resulting $Det(\overline{C})=1$, constant.

Defining

$$\widehat{C}={\widehat{A}}^T\widehat{A}=J^{2/3}I$$,(20)

results in the desired multiplicative decomposition applied to the Cauchy-Green stretch tensor as:

$$C=\widehat{C}\cdot \overline{C}=\overline{C}\cdot \widehat{C}$$ (21)

One may observe that $\det (\widehat{C})=J^2=\det (C)$. From $\overline{C}$ and $\widehat{C}$ one writes and improves various hyperelastic constitutive models, see for instance Rivlin and Saunders (1951), Bonet and Wood (1997), Düster et al. (2003), Hartmann and Neff (2003).

Lagrangian finite strain directions definition

The definition and understanding of the Lagrangian finite strain directions (volumetric and isochorics) start by the analysis of the derivative of parts of the hyperelastic Helmholtz free energy (Eq. (15)) regarding the Green-Lagrange strain.

Observing Equation (20) and comments after Equation (21), the volumetric elastic part of the Helmholtz free energy potential (Eq. (15) is written as an exclusive function of the Jacobian of the deformation function, as:

$${\psi }_{elas}^{vol}={\psi }_{elas}^{vol}\left(J\right)={\psi }_{elas}^{vol}\left(\sqrt{\det (\widehat{C})}\right)$$ (22)

from which one calculates the corresponding second Piola-Kirchhoff stress as:

$$S_{elas}^{vol}=\frac{\partial {\psi }_{elas}^{vol}}{\partial E}=\frac{\partial {\psi }_{elas}^{vol}}{\partial J}\frac{\partial J}{\partial E}=\alpha {\mathfrak{E}}^{vol}$$ (23)

with $\alpha =\partial {\psi }_{elas}^{vol}/\partial J$ being a scalar and ${\mathfrak{E}}^{vol}$ a symmetric tensor of the same order of the Green strain. It is simple to show that:

$${\mathfrak{E}}^{vol}=\frac{\partial J}{\partial E}=J{C}^{-1}$$ (24)

Applying the first part of Eq. (12) over the second Piola-Kirchhoff stress of Eq. (23) (pulling forward) results in:

$${\sigma }_{elas}^{vol}=\frac{1}{J}A{S}_{elas}^{vol}{A}^T=\alpha J\frac{1}{J}A{C}^{-1}{A}^T=\alpha I={\sigma }_{hyd}$$,(25)

which means that ${\mathfrak{E}}_{vol}$ is the Lagrangian strain direction corresponding to the hydrostatic stress in the Cauchy space and $\alpha =\partial {\psi }_{elas}^{vol}/\partial J$ corresponds to the hydrostatic second Piola-Kirchhoff stress intensity. The identification of the Lagrangian hydrostatic stress $S_{elas}^{vol}$ and the corresponding strain direction ${\mathfrak{E}}^{vol}$ is important to define the proposed viscous model.

The first isochoric elastic part of the Helmholtz free energy potential is written as an exclusive function of the first invariant (${\overline{I}}_1$) of the isochoric part of the Cauchy-Green stretch tensor $\overline{C}$ as:

$${\psi }_{elas}^{iso1}={\psi }_{elas}^{iso1}\left({\overline{I}}_1\right)$$,(26)

from which one calculates the corresponding second Piola-Kirchhoff stress as:

$$S_{elas}^{iso1}=\frac{\partial {\psi }_{elas}^{iso1}}{\partial {\overline{I}}_1}\frac{\partial {\overline{I}}_1}{\partial E}=\beta {\mathfrak{E}}^{iso1}$$ (27)

where $\beta$ is a scalar and ${\mathfrak{E}}^{iso1}$ is a symmetric tensor of the same order of the Green strain. Opening ${\overline{I}}_1$ one writes:

$${\mathfrak{E}}^{iso1}=\frac{\partial {\overline{I}}_1}{\partial E}=\frac{\partial tr\left(\overline{C}\right)}{\partial E}=\frac{\partial \left(J^{-2/3}tr\left(C\right)\right)}{\partial E}$$ (28)

from which one finds:

$${\mathfrak{E}}^{iso1}=-\frac{2}{3}{J}^{-2/3}tr(C)C^{-1}+2J^{-2/3}I$$ (29)

Considering Eq. (29) and applying the transformation (12) over Eq. (27) one finds:

$${\sigma }_{iso1}=\beta \left\{2J^{-5/3}\left(A{A}^t-\frac{tr\left(A^{T}A\right)}{3}I\right)\right\}$$ (30)

As $$tr\left(A^{T}A\right)=A:A^T=A^T:A=tr(A{A}^T)$$, one concludes that $${\sigma }_{iso}$$ is a deviatory $${\sigma }_{dev}$$ stress, i.e.:

$${\sigma }_{iso1}=\beta \left\{2J^{-5/3}\left(A{A}^T-\frac{tr\left(A{A}^T\right)}{3}I\right)\right\}={\sigma }_{dev}$$ (31)

which means that ${\mathfrak{E}}^{iso1}/\sqrt{{\mathfrak{E}}^{iso1}:{\mathfrak{E}}^{iso1}}$ is the first isochoric Lagrangian strain direction corresponding to a deviatory strain, and $\beta =\left(\partial {\psi }_{elas}^{iso1}/\partial {\overline{I}}_1\right)\sqrt{{\mathfrak{E}}^{iso1}:{\mathfrak{E}}^{iso1}}$ is the associated second Piola-Kirchhoff stress intensity.

The second isochoric elastic part of the Helmholtz free energy potential is written as an exclusive function of the second invariant of the isochoric Cauchy-Green stretch tensor ${\overline{I}}_2$ as:

$${\psi }_{elas}^{iso2}={\psi }_{elas}^{iso2}\left({\overline{I}}_2\right)$$ (32)

from which we calculate the corresponding second Piola-Kirchhoff stress as:

$$S_{elas}^{iso2}=\frac{\partial {\psi }_{elas}^{iso2}}{\partial {\overline{I}}_2}\frac{\partial {\overline{I}}_2}{\partial E}=\gamma {\mathfrak{E}}^{iso2}$$ (33)

with $\gamma$ and ${\mathfrak{E}}^{iso2}$being respectively a scalar and a symmetric tensor of the same order as the Green-Lagrange strain tensor. It is interesting to write:

$${\overline{I}}_2=J^{-4/3}{I}_2=J^{-4/3}\left\{\left(C_{11}{C}_{22}-C_{12}{C}_{21}\right)+\left(C_{11}{C}_{33}-C_{13}{C}_{31}\right)+\left(C_{22}{C}_{33}-C_{23}{C}_{32}\right)\right\}$$ (34)

in which $I_2$ is the second invariant of the Cauchy stretch. With some algebraic effort over Eq. (34), results:

$${\mathfrak{E}}^{iso2}=\frac{\partial {\overline{I}}_2}{\partial E}=2J^{-4/3}\left(-\frac{2}{3}{C}^{-1}{I}_2+\left\{tr(C)I-C^T\right\}\right)$$ (35)

Using Eq. (12) over Eq. (33) and considering Eq. (35), one writes:

$${\sigma }_{iso2}=\gamma 2J^{-7/3}\left[\left(tr(C)\left(A{A}^T\right)-\left(A{A}^T\right)\left(A{A}^T\right)\right)-\left(\frac{2}{3}{I}_2\right)I\right]$$ (36)

Again, as $tr\left(C\right)=tr\left(A^{T}A\right)=tr\left(A{A}^T\right)$, results:

$$tr(A{A}^T)tr(A{A}^T)-tr((A{A}^T)(A{A}^T))=2I_2$$ (37)

From Equation (36), using Eq. (37), one calculates,

$$tr\left({\sigma }_{iso2}\right)=2\gamma J^{-7/3}\left(2I_2-2I_2\right)=0$$ (38)

thus:

$${\sigma }_{iso2}={\sigma }_{dev}$$ (39)

In equation (39) one concludes that ${\sigma }_{iso2}$ is deviatory and, thus ${\mathfrak{E}}^{iso2}/\sqrt{{\mathfrak{E}}^{iso2}:{\mathfrak{E}}^{iso2}}$ is the second isochoric Lagrangian strain direction, and $\gamma =\left(\partial {\psi }_{elas}^{iso2}/\partial {\overline{I}}_2\right)\sqrt{{\mathfrak{E}}^{iso2}:{\mathfrak{E}}^{iso2}}$ is the intensity of the second deviatory Piola-Kirchhoff stress. Thus, ${\mathfrak{E}}_{vol}$, ${\mathfrak{E}}^{iso1}/\sqrt{{\mathfrak{E}}^{iso1}:{\mathfrak{E}}^{iso1}}$ and ${\mathfrak{E}}^{iso2}/\sqrt{{\mathfrak{E}}^{iso2}:{\mathfrak{E}}^{iso2}}$ are the Lagrangian version of hydrostatic, deviatory 1 and deviatory 2 uncoupled finite strain directions.

The adopted Hyperelastic part of the model

We chose a Rivlin-Saunders-Düster-Hartmann type hyperelastic constitutive model, Rivlin and Saunders (1951) and Düster et al. (2003), i.e.:

$${\psi }_{vol}=\frac{K}{8n^2}\left(J^{2n}+J^{-2n}-2\right)$$ (40)

and

$${\psi }_{iso1}^{}+{\psi }_{iso2}^{}=\frac{G}{4}\left({\overline{I}}_1-3\right)+\frac{G}{4}\left({\overline{I}}_2-3\right)$$ (41)

In which $K$ is the bulk modulus, $G$ is the transversal elastic modulus and $n\ne 0$ is an integer. Other hyperelastic potentials that respect the volumetric and isochoric split can also be chosen.

Differentiating equations (40) and (41) with respect to Green strain, using information and results of equations (15), (23), (27) and (33), we find the second Piola-Kirchhoff elastic stress as:

$$S_{elas}=\left\{\frac{K}{4}\left(J-J^{-3}\right)\right\}{\mathfrak{E}}^{vol}+\frac{G}{4}{\mathfrak{E}}^{iso1}+\frac{G}{4}{\mathfrak{E}}^{iso2}$$ (42)

in which $n=1$ is adopted here.

2.3 Proposed simple Kelvin-Voigt-like solids constitutive model including viscous fluids Lagrangian modeling

Let us consider the Helmholtz free energy (scalar) written as the sum of two parts, one hyperelastic (Equation(42)) and another viscous. This sum can only be represented in its differential form (Lánczos (1970) and Sanches and Coda HB (2013)), as there is not an explicit representation for dissipation, i.e.:

$$\delta \psi =\delta {\psi }_{elas}+\delta {\psi }_{vis}=S_{elas}^{vol}:\delta E+S_{elas}^{iso1}:\delta E+S_{elas}^{iso2}:\delta E+S_{vis}^{vol}:\delta E+S_{vis}^{iso1}:\delta E+S_{vis}^{iso2}:\delta E$$ (43)

or, splitting

$$S:\delta E=S_{elas}^{}:\delta E+S_{vis}^{}:\delta E=\left(\frac{\partial {\psi }_{elas}^{vol}}{\partial E}+\frac{\partial {\psi }_{elas}^{iso1}}{\partial E}+\frac{\partial {\psi }_{elas}^{iso2}}{\partial E}\right):\delta E+\left(S_{vis}^{vol}+S_{vis}^{iso1}+S_{vis}^{iso2}\right):\delta E$$ (44)

in which the elastic stress components are given by Eq. (42) from the definition given in the second of Eqs. (14) and the steps described by Eqs. (23), (27) and (33).

The primary idea to write the viscous stress is to make an analogous of the hyperelastic relation (Eq. (42)) to incorporate viscosity assuming that the time rate of the split strain directions would suite a constitutive relation, i.e.:

$$\delta {\psi }_{vis}=S_{vis}^{*}:\delta E=\left(\frac{\overline{K}}{4}{\dot{\mathfrak{E}}}^{vol}+\frac{\overline{G}}{4}{\dot{\mathfrak{E}}}^{iso1}+\frac{\overline{G}}{4}{\dot{\mathfrak{E}}}^{iso2}\right):\delta E$$ (45)

However, contrary to small strains (Mesquita and Coda (2007a) and Mesquita and Coda (2007b)), the time derivatives of volumetric and isochoric finite strain directions do not preserve direction. Thus, $S_{vis}^{*}$ serves only as an inspiration to the following developments.

Inspired in Eqs. (40,41,42,45)))))))) in order to keep isotropy, the viscous virtual work should be written regarding the proper scalar variables with their variations expanded to comprise strains, i.e., $\delta J=\partial J/\partial E:\delta E$ and $\delta {\overline{I}}_i=\partial {\overline{I}}_i/\partial E:\delta E$. This is done by:

$$\begin{array}{l}\delta \psi =\frac{\overline{K}}{4\sqrt{{\mathfrak{E}}^{vol}:{\mathfrak{E}}^{vol}}}\frac{d{J}^{\alpha }}{dt}\delta J+\frac{{\overline{G}}_{(i)}}{4\sqrt{{\mathfrak{E}}^{iso(i)}:{\mathfrak{E}}^{iso(i)}}}\frac{d{\overline{I}}_i^{\gamma }}{dt}\delta {\overline{I}}_i=\\ =\left(\frac{\overline{K}}{4\sqrt{{\mathfrak{E}}^{vol}:{\mathfrak{E}}^{vol}}}\alpha J^{\alpha -1}\dot{J}\frac{\partial J}{\partial E}+\frac{{\overline{G}}_1}{4\sqrt{{\mathfrak{E}}^{iso1}:{\mathfrak{E}}^{iso1}}}{\gamma }_1{\overline{I}}_1^{{\gamma }_1-1}{\dot{\overline{I}}}_1\frac{\partial {\overline{I}}_1}{\partial E}+\frac{{\overline{G}}_2}{4\sqrt{{\mathfrak{E}}^{iso2}:{\mathfrak{E}}^{iso2}}}{\gamma }_2{\overline{I}}_2^{{\gamma }_2-1}{\dot{\overline{I}}}_2\frac{\partial {\overline{I}}_2}{\partial E}\right):\delta E\end{array}$$ (46)

From Eq. (46) and Eqs. (24), (29) and (35) the following expression for the second Piola-Kirchhoff viscous stress is written,

$$S_{vis}=\frac{\overline{K}}{4}\frac{\alpha J^{\alpha -1}\dot{J}{\mathfrak{E}}^{vol}}{\sqrt{{\mathfrak{E}}^{vol}:{\mathfrak{E}}^{vol}}}+\frac{{\overline{G}}_1}{4}\frac{{\gamma }_1{\overline{I}}_1^{{\gamma }_1-1}{\dot{\overline{I}}}_1{\mathfrak{E}}^{iso1}}{\sqrt{{\mathfrak{E}}^{iso1}:{\mathfrak{E}}^{iso1}}}+\frac{{\overline{G}}_2}{4}\frac{{\gamma }_2{\overline{I}}_2^{{\gamma }_2-1}{\dot{\overline{I}}}_2{\mathfrak{E}}^{iso2}}{\sqrt{{\mathfrak{E}}^{iso2}:{\mathfrak{E}}^{iso2}}}$$ (47)

in which $\overline{K}$ is the fluid volumetric viscosity Holmes et al. (2011) and ${\overline{G}}_i$ are shear (isochoric directions 1 and 2) viscosities and ${\mathfrak{E}}^{vol}$, ${\mathfrak{E}}^{iso1}$ and ${\mathfrak{E}}^{iso2}$ are strain directions. In order to be coherent with the elasticity understanding we assume $\overline{G}={\overline{G}}_1={\overline{G}}_2$. Adopting the viscous parameters ${\gamma }_1={\gamma }_2=1/2$ and $\overline{K}=0$ one reproduces Newtonian fluids and a simple Kelvin-Voigt visco-hyperelasticity. Notice that, as time strain rates are written as function of dimensionless scalars (strain invariants), the values of $\alpha$ and ${\gamma }_i$ can be investigated for different viscous behavior. When these parameters are null, logarithm viscosity is assumed. The investigation of the fluid behavior regarding different choices of $\alpha$ and ${\gamma }_i$ is out of this study objectives; however some inspections are made in example 4.5.

In short, Eq. (44) becomes:

$$S:\delta E=\left(S_{elas}+S_{vis}\right):\delta E$$ or $$S=S_{elas}+S_{vis}$$ (48)

For fluids, when calculating $S_{elas}$ only the volumetric part of Eq. (42) is considered, i.e. $G=0$. When modeling elastic solids $G\ne 0$ and $K\ne 0$ but $\overline{G}=\overline{K}=0$. For viscoelastic solids all elastic and viscous constants are taken into account. Moreover, it is important to note that the volumetric viscosity is disregarded $\overline{K}=0$ unless for specific cases Holmes et al. (2011).

To define the numerical solution procedure (in the next item) it is necessary to know the viscoelastic tangent constitutive tensor, given by:

$$ℭ=\frac{\partial S}{\partial E}=\frac{\partial S_{elas}}{\partial E}+\frac{\partial S_{vis}}{\partial E}={ℭ}_{elas}+{\mathfrak{N}}_{vis}$$ (49)

The first term can be calculated as:

$${ℭ}_{elas}={ℭ}^{vol}+{ℭ}^{iso1}+{ℭ}^{iso2}=\frac{{\partial }^2{\Psi }_{elas}}{\partial E\otimes \partial E}$$ (50)

with

$${ℭ}_{}^{vol}=\frac{{\partial }^2{\Psi }_{vol}}{\partial E\otimes \partial E}=\frac{{\partial }^2{\Psi }_{vol}}{\partial J^2}\frac{\partial J}{\partial E}\otimes \frac{\partial J}{\partial E}+\frac{\partial {\Psi }_{vol}}{\partial J}\frac{{\partial }^{2}J}{\partial E\otimes \partial E}$$ (51)

$${ℭ}^{iso1}=\frac{{\partial }^2{\Psi }_{iso1}}{\partial E\otimes \partial E}=\frac{{\partial }^2{\Psi }_{iso}^1}{{\left(\partial {\overline{I}}_1\right)}^2}\frac{\partial {\overline{I}}_1}{\partial E}\otimes \frac{\partial {\overline{I}}_1}{\partial E}+\frac{\partial {\Psi }_{iso}^1}{\partial {\overline{I}}_1}\frac{{\partial }^2{\overline{I}}_1}{\partial E\otimes \partial E}$$ (52)

$${ℭ}^{iso2}=\frac{{\partial }^2{\Psi }_{iso2}}{\partial E\otimes \partial E}=\frac{{\partial }^2{\Psi }_{iso}^2}{{\left(\partial {\overline{I}}_2\right)}^2}\frac{\partial {\overline{I}}_2}{\partial E}\otimes \frac{\partial {\overline{I}}_2}{\partial E}+\frac{\partial {\Psi }_{iso}^2}{\partial {\overline{I}}_2}\frac{{\partial }^2{\overline{I}}_2}{\partial E\otimes \partial E}$$ (53)

in which $\partial J/\partial E$ is given by Eq. (24), $\partial {\overline{I}}_1/\partial E$ is given by Eq. (29), $\partial {\overline{I}}_2/\partial E$ is given by Eq. (35) and the other terms are given, as:

$$\frac{{\partial }^2{\Psi }_{vol}}{\partial J^2}=\frac{K}{4}\left(1+3J^{-4}\right)$$ (54)

$$\frac{{\partial }^{2}J}{\partial E_{j\gamma }\partial E_{oz}}=J\left\{D_{j\gamma }{D}_{oz}-2D_{jo}{D}_{z\gamma }\right\}=\frac{\partial {\mathfrak{E}}_{j\gamma }^{vol}}{\partial E_{oz}}$$ (55)

$$\frac{{\partial }^2{\overline{I}}_1}{\partial E_{ij}\partial E_{k\mathcal{l}}}=\frac{4}{3}{J}^{-2/3}\left\{\frac{1}{3}\left\{D_{ij}{D}_{k\mathcal{l}}+3D_{ik}{D}_{\mathcal{l}j}\right\}{I}_1-D_{ij}{\delta }_{k\mathcal{l}}-D_{k\mathcal{l}}{\delta }_{ij}\right\}=\frac{\partial {\mathfrak{E}}_{ij}^{iso1}}{\partial E_{k\mathcal{l}}}$$ (56)

and

$$\frac{{\partial }^2{\overline{I}}_2}{\partial E_{ij}\partial E_{k\mathcal{l}}}=\frac{8}{3}{J}^{-4/3}\left\{\begin{array}{l}\left[\frac{2}{3}{D}_{ij}{D}_{k\mathcal{l}}+D_{ik}{D}_{\mathcal{l}j}\right]I_2-\left[C_{zz}\left(D_{ij}{\delta }_{k\mathcal{l}}+D_{k\mathcal{l}}{\delta }_{ij}\right)\right]+\\ +D_{ij}{C}_{k\mathcal{l}}+D_{k\mathcal{l}}{C}_{ij}+\frac{3}{2}\left[{\delta }_{ij}{\delta }_{k\mathcal{l}}-{\delta }_{jk}{\delta }_{i\mathcal{l}}\right]\end{array}\right\}=\frac{\partial {\mathfrak{E}}_{ij}^{iso2}}{\partial E_{k\mathcal{l}}}$$ (57)

in which $D_{ij}=C_{ij}^{-1}$. One may note that in this model ${\partial }^2{\Psi }_{iso}^1/{\left(\partial {\overline{I}}_1\right)}^2=0$ and ${\partial }^2{\Psi }_{iso}^2/{\left(\partial {\overline{I}}_2\right)}^2=0$.

As the goal of this study is to give a numerical solution for both fluids and solids, we approximate the strain invariant rates in Eq. (47) by backward finite difference (usual in fluid mechanics), instead of using backward-Euler stress integration (usual in solid mechanics) so that:

$$S_{vis}=\frac{\overline{K}}{4}\alpha J_{s+1}^{\alpha -1}\left(\frac{J_{s+1}-J_s}{\Delta t}\right){\mathfrak{E}}^{vol}+\frac{{\overline{G}}_1}{4}{\gamma }_1{\overline{I}}_{1(s+1)}^{{\gamma }_1-1}\left(\frac{{\left({\overline{I}}_1\right)}_{s+1}-{\left({\overline{I}}_1\right)}_s}{\Delta t}\right){\mathfrak{E}}^{iso1}+\frac{{\overline{G}}_2}{4}{\gamma }_2{\overline{I}}_{2(s+1)}^{{\gamma }_2-1}\left(\frac{{\left({\overline{I}}_2\right)}_{s+1}-{\left({\overline{I}}_2\right)}_s}{\Delta t}\right){\mathfrak{E}}^{iso2}$$ (58)

Differentiating Eq. (58) regarding current strain (implicit method), it results the following numeric tangent constitutive tensor at the current time ($s+1$):

$$\begin{array}{l}{\mathfrak{N}}_{}^{vis}=\frac{\partial S_{vis}}{\partial E}=\frac{\overline{K}}{4\Delta t}\left({\alpha }^2{J}_{}^{\alpha -1}{\mathfrak{E}}^{vol}\otimes {\mathfrak{E}}^{vol}+\alpha J_{}^{\alpha }\frac{{\partial }^{2}J}{\partial E\otimes \partial E}\right)+\\ +\frac{{\overline{G}}_1}{4\Delta t}\left({\gamma }_1^2{\overline{I}}_1^{{\gamma }_1-1}{\mathfrak{E}}^{iso1}\otimes {\mathfrak{E}}^{iso1}+{\gamma }_1{\overline{I}}_1^{{\gamma }_1}\frac{{\partial }^2{\overline{I}}_1^{}}{\partial E\otimes \partial E}\right)+\\ +\frac{{\overline{G}}_2}{4\Delta t}\left({\gamma }_2^2{\overline{I}}_2^{{\gamma }_2-1}{\mathfrak{E}}^{iso2}\otimes {\mathfrak{E}}^{iso2}+{\gamma }_2{\overline{I}}_2^{{\gamma }_2}\frac{{\partial }^2{\overline{I}}_2^{}}{\partial E\otimes \partial E}\right)\end{array}$$ (59)

Notice that terms related to the past (of Eq. (58)) do not appear in Eq. (59) as they are constant. Substituting ${\mathfrak{E}}^{vol}$, ${\mathfrak{E}}^{iso1}$ and ${\mathfrak{E}}^{iso2}$ by their respective definitions (24), (29) and (35), and using index notation, it results:

$$\begin{array}{l}{\mathfrak{N}}_{ijk\mathcal{l}}^{vis}=\frac{1}{4\Delta t}\left\{\overline{K}\left({\alpha }^2{J}_{}^{\alpha -1}\frac{\partial J}{\partial E_{ij}}\frac{\partial J}{\partial E_{k\mathcal{l}}}+\alpha J_{}^{\alpha }\frac{{\partial }^{2}J}{\partial E_{ij}\partial E_{k\mathcal{l}}}\right)\right\}+\\ +{\overline{G}}_1\left({\gamma }_1^2{\overline{I}}_1^{{\gamma }_1-1}\frac{\partial {\overline{I}}_1^{}}{\partial E_{ij}}\frac{\partial {\overline{I}}_1^{}}{\partial E_{k\mathcal{l}}}+{\gamma }_1{\overline{I}}_1^{{\gamma }_1}\frac{{\partial }^2{\overline{I}}_1^{}}{\partial E_{ij}\partial E_{k\mathcal{l}}}\right)+\\ \left.+{\overline{G}}_2\left({\gamma }_2^2{\overline{I}}_2^{{\gamma }_1-1}\frac{\partial {\overline{I}}_2^{}}{\partial E_{ij}}\frac{\partial {\overline{I}}_2^{}}{\partial E_{k\mathcal{l}}}+{\gamma }_2{\overline{I}}_2^{{\gamma }_2}\frac{{\partial }^2{\overline{I}}_2^{}}{\partial E_{ij}\partial E_{k\mathcal{l}}}\right)\right\}\end{array}$$ (60)

in which, for simplicity, the current time is not indicated.

3 BRIEF DESCRIPTION OF THE IMPLEMENTED POSITIONAL FEM

Position-based FEM has been used by the authors in several applications, see Coda (2015), Coda HB and Paccola (2011) and Pascon and Coda (2013) for instance. The difference between this alternative and the classic FEM is the adoption of positions as nodal parameters instead of displacements making solid mechanics description more compact and direct. This gives to the technique a simple association with geometric nonlinear problems of deformable solids that is easily extended here to model low-compressible free-surface flows.

Firstly, one chooses the nodal positions $\overrightarrow{Y}$ as the analysis parameters, thus the continuum position and acceleration are represented by:

$$\overrightarrow{y}=\varphi \cdot \overrightarrow{Y}$$ (61)

$$\overrightarrow{\ddot{y}}=\varphi \cdot \overrightarrow{\ddot{Y}}$$ (62)

in which $\varphi$ are usual FEM shape functions. Volume and surface forces, if desired, can be approximated by their nodal values, as:

$${\overrightarrow{b}}^0=\varphi \cdot {\overrightarrow{B}}^0$$ (63)

and

$$\overrightarrow{h}=\phi \cdot \overrightarrow{P}$$ (64)

where $\phi$ are surface shape functions.

The strain variation $\delta E$ should also be given as a function of position variation $\delta \overrightarrow{Y}$, as:

$$\delta E=\frac{\partial E}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}$$ (65)

Introducing such changes and considereing Eq. (48), the weak form of the equilibrium Eq. (13) becomes:

$$\left({\displaystyle {\int }_{V_0}^{}{\rho }_0\varphi \otimes \varphi d{V}_0\cdot \overrightarrow{\ddot{Y}}}-{\displaystyle {\int }_{V_0}^{}\varphi \otimes \varphi d{V}_0\cdot {\overrightarrow{B}}^0}-{\displaystyle {\int }_{S_0}^{}\phi \otimes \phi d{S}_0\cdot \overrightarrow{P}}+{\displaystyle {\int }_{V_0}^{}\left(S_{elas}+S_{vis}\right):\frac{\partial E}{\partial \overrightarrow{Y}}d{V}_0}\right)\cdot \delta \overrightarrow{Y}=0$$ (66)

As the nodal position variation $\delta \overrightarrow{Y}$ is arbitrary, Eq. (66) turns into the following nonlinear system of equations:

$$M\cdot \overrightarrow{\ddot{Y}}-{\overrightarrow{F}}^{ext}+{\overrightarrow{F}}^{int}=\overrightarrow{0}$$ (67)

in which the last term of Eq. (66) has been called internal force and its (nonlinear) calculation is clarified in the definition of the finite element kinematics. Concentrated external loads can be applied directly in ${\overrightarrow{F}}^{ext}$.

In this study, the continuum (fluids and solids) is discretized by triangular-base prismatic three-dimensional finite elements with linear or cubic approximations in their base directions (${\xi }_1,{\xi }_2$), and any order of approximation along its thickness direction (${\xi }_3$). Fig. 1 shows a finite element with linear approximation at ${\xi }_3$ direction and cubic approximation at ${\xi }_1$ and ${\xi }_2$ directions.

Fig. 1
Deformation function described by prismatic elements

The position based kinematics description is given in a brief way as more details can be seen in references Sanches and Coda (2013), Coda (2015), Coda HB and Paccola (2011) and Pascon and Coda (2013) for instance. In Fig. 1 one observes that the initial configuration mapping describes the initial coordinates of continuous points $x_i$ as a function of the dimensionless coordinates ${\xi }_i$ and of the initial coordinates $X_{\alpha i}$ of nodes $\alpha$. Similarly, current continuum coordinates $y_i$ are written as function of the same dimensionless coordinates and current nodal coordinates (unknown of the problem) $Y_{\alpha i}$. These mappings are represented in index notation as:

$$f_i^0=x_i={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right)X_{\alpha i}$$ and $$f_i^1=y_i={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right)Y_{\alpha i}$$ (68)

in which the last one is an index version of Eq. (61). The gradients of these mappings are given by:

$$A_{ij}^0=\frac{\partial x_i}{\partial {\xi }_j}$$ and $$A_{ij}^1=\frac{\partial y_i}{\partial {\xi }_j}$$ (69)

Thus, the gradient $A$ of the deformation function $\overrightarrow{f}$ (see Equations (12) and (16)) is given by:

$$A=A^1{\left(A^0\right)}^{-1}$$ (70)

Eq. (70) is numerically written at each integration point for a trial position $\overrightarrow{Y}$ covering all nodes of the problem. As the formulation is total Lagrangian, one may observe that the gradient of the initial configuration mapping $A^0$ is calculated only once in the numerical process. Moreover, it is clear the dependence of $A$ regarding $\overrightarrow{Y}$.

After calculating the deformation gradient $A$, one calculates the Cauchy-Green stretch tensor and the Green strain as:

$$C=A^{T}A$$ and $$E=\frac{1}{2}\left(C-I\right)$$ (71)

Using the equations of the previous sections, all necessary variables, such as ${\mathfrak{E}}^{vol}$, ${\mathfrak{E}}^{iso1}$, ${\mathfrak{E}}^{iso2}$, $\dot{J}$, ${\dot{\overline{I}}}_1$ and ${\dot{\overline{I}}}_2$, can be numerically calculated at each integration point of Eq. (66) or Eq. (67) and are function of nodal positions $\overrightarrow{Y}$.

Eq. (67) is valid at any instant and, in the numerical process; the current instant is written as the previous instant plus a time step, as:

$$t_{s+1}=t_s+\Delta t$$ (72)

in which $t_{s+1}$ is the current instant. Thus, one writes Eq. (67) as:

$${\overrightarrow{g}}_{s+1}=M\cdot {\overrightarrow{\ddot{Y}}}_{s+1}+C\cdot {\overrightarrow{\dot{Y}}}_{s+1}-{\overrightarrow{F}}_{s+1}^{ext}+\left({\overrightarrow{F}}_{s+1}^{vis}+{\overrightarrow{F}}_{s+1}^{elas}\right)=\overrightarrow{0}$$ (73)

in which, for simplicity, time $t$ is omitted and the internal force is split in elastic and viscous parts. For completeness, a damping matrix $C$ is also included; following usual solid mechanic works Coda (2015).

As the formulation is total Lagrangian (constant mass matrix), we adopted the Newmark approximation to perform the time marching, written as:

$${\overrightarrow{Y}}_{s+1}={\overrightarrow{Y}}_s+{\overrightarrow{\dot{Y}}}_s\Delta t+\left[\left(\frac{1}{2}-\beta \right){\overrightarrow{\ddot{Y}}}_s+\beta {\overrightarrow{\ddot{Y}}}_{s+1}\right]\Delta t^2$$ (74)

$${\overrightarrow{\dot{Y}}}_{s+1}={\overrightarrow{\dot{Y}}}_s+\left(1-\gamma \right)\Delta t{\overrightarrow{\ddot{Y}}}_s+\gamma \Delta t{\overrightarrow{\ddot{Y}}}_{s+1}$$ (75)

in which $\beta$ and $\gamma$ are free parameters, generally adopted as $\beta =1/4$ and $\gamma =1/2$.

Isolating the current velocity and the current acceleration from Eqs. (74) and (75), results:

$${\overrightarrow{\ddot{Y}}}_{s+1}=\frac{{\overrightarrow{Y}}_{s+1}}{\beta \Delta t^2}-{\overrightarrow{Q}}_s$$ (76)

and

$${\dot{\overrightarrow{Y}}}_{s+1}=\frac{\gamma }{\beta \Delta t}{\overrightarrow{Y}}_{s+1}+{\overrightarrow{R}}_s-\gamma \Delta t{\overrightarrow{Q}}_s$$ (77)

with the following auxiliary values:

$${\overrightarrow{Q}}_s=\left(\frac{{\overrightarrow{Y}}_s}{\beta \Delta t^2}+\frac{{\dot{\overrightarrow{Y}}}_s}{\beta \Delta t^{}}+\left(\frac{1}{2\beta }-1\right){\ddot{Y}}_s\right)$$ and $${\overrightarrow{R}}_s=\left[{\dot{\overrightarrow{Y}}}_s+\Delta t\left(1-\gamma \right){\ddot{\overrightarrow{Y}}}_s\right]$$ (78)

Substituting Eqs. (76) and (77) into (73), results:

$$\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)={\overrightarrow{F}}_{s+1}^{elas}+{\overrightarrow{F}}_{s+1}^{vis}+\frac{M}{\beta \Delta t^2}{\overrightarrow{Y}}_{s+1}-M\cdot {\overrightarrow{Q}}_s\text{+}\frac{\gamma C}{\beta \Delta t}\cdot {\overrightarrow{Y}}_{s+1}+C\cdot {\overrightarrow{R}}_s-\gamma \Delta tC\cdot {\overrightarrow{Q}}_s-{\overrightarrow{F}}_{s+1}^{ext}(t)=\overrightarrow{0}$$ (79)

Eq. (79) is understood simply as $\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)=0$, revealing the nonlinear behavior of equilibrium equations regarding $\left({\overrightarrow{Y}}_{s+1}\right)$, solved in this work by the Newton Raphson method.

Following previous works Coda (2015), Coda HB and Paccola (2011) and Pascon and Coda (2013), using a truncated Taylor expansion, we have:

$$\overrightarrow{0}=\overrightarrow{g}({\overrightarrow{Y}}_{s+1}^{})\cong \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)+\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)\Delta \overrightarrow{Y}$$ (80)

where

$$\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}\right)=H={\left.\frac{\partial {\overrightarrow{F}}^{elas}}{\partial \overrightarrow{Y}}\right|}_{s+1}+{\left.\frac{\partial {\overrightarrow{F}}^{vis}}{\partial \overrightarrow{Y}}\right|}_{s+1}+\left(\frac{M}{\beta \Delta t^2}+\frac{\gamma C}{\beta \Delta t}\right)=H^{elas}+H^{vis}+H^{dyn}$$ (81)

with $H^{elas}$ and $H^{din}$ being, respectively, the elastic and dynamic parts of the Hessian matrix. The viscous term $H^{vis}$ is an original expression introduced in this study and is a key value in the total Lagrangian isothermal fluid model and simple viscoelastic model proposed here, see Eq. (90).

From Eq. (80), we write the linear system of equation from which the position correction $\Delta \overrightarrow{Y}$ is achieved as:

$$\nabla \overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)\cdot \Delta \overrightarrow{Y}=-\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)$$ or $$H\cdot \Delta \overrightarrow{Y}=-\overrightarrow{g}\left({\overrightarrow{Y}}_{s+1}^0\right)$$ (82)

in which ${\overrightarrow{Y}}_{s+1}^0$ is a trial position. In the beginning of a time step ${\overrightarrow{Y}}_{s+1}^0$ is assumed as the result of the previous step, i.e., ${\overrightarrow{Y}}_s^{}$ and at the end of the time step ${\overrightarrow{Y}}_{s+1}^0$ becomes the accepted solution ${\overrightarrow{Y}}_{s+1}^{}$. Solving the correction $\Delta \overrightarrow{Y}$ in Eq. (82), a new trial for ${\overrightarrow{Y}}_{s+1}$ is calculated as:

$${\overrightarrow{Y}}_{s+1}^0={\overrightarrow{Y}}_{s+1}^0+\Delta \overrightarrow{Y}$$ (83)

At the beginning of each iteration, acceleration and velocity should be recalculated as:

$${\overrightarrow{\ddot{Y}}}_{s+1}=\frac{{\overrightarrow{Y}}_{s+1}^0}{\beta \Delta t^2}-{\overrightarrow{Q}}_S$$ and $${\dot{\overrightarrow{Y}}}_{s+1}=\frac{\gamma }{\beta \Delta t}{\overrightarrow{Y}}_{s+1}^0+{\overrightarrow{R}}_s-\gamma \Delta t{\overrightarrow{Q}}_s$$ (84)

As ${\overrightarrow{Q}}_s$ and ${\overrightarrow{R}}_s$ are values of the past, they remain unaltered during iterations. The stopping test is given by:

$$\frac{‖\overrightarrow{g}({\overrightarrow{Y}}^0)‖}{‖{\overrightarrow{F}}^{ext}‖}\le TOL$$ or $$\frac{‖\Delta \overrightarrow{Y}‖}{‖\overrightarrow{X}‖}\le TOL$$ (85)

in which $TOL$ is the prescribed numerical tolerance.

At the first time step, the initial acceleration is calculated by the dynamic equilibrium equation as:

$${\overrightarrow{\ddot{Y}}}_0=M^{-1}\left[{\overrightarrow{F}}_0^{ext}-{\overrightarrow{F}}^{elas}-{\overrightarrow{F}}^{vis}-C\cdot {\overrightarrow{\dot{Y}}}_0\right]$$ (86)

The elastic and viscous parts of Hessian matrix are derived for each finite element from the derivative of Eq. (73) (see Eq. (81)) regarding the current position vector, as:

$$H^{elas}=\frac{\partial {\overrightarrow{F}}^{elas}}{\partial \overrightarrow{Y}}={\displaystyle {\int }_{V_0}^{}\frac{\partial }{\partial \overrightarrow{Y}}\left(\frac{\partial {\psi }_{elas}}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}\right)d{V}_0}={\displaystyle {\int }_{V_0}^{}\left(\frac{\partial E}{\partial \overrightarrow{Y}}:\frac{{\partial }^2{\psi }_{elas}}{\partial E\otimes \partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}+\frac{\partial {\psi }_{elas}}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\right)d{V}_0}$$ (87)

$$H^{vis}=\frac{\partial {\overrightarrow{F}}^{vis}}{\partial \overrightarrow{Y}}={\displaystyle {\int }_{V_0}^{}\frac{\partial }{\partial \overrightarrow{Y}}\left(S_{vis}:\frac{\partial E}{\partial \overrightarrow{Y}}\right)d{V}_0}={\displaystyle {\int }_{V_0}^{}\left(\frac{\partial E}{\partial \overrightarrow{Y}}:\frac{\partial S_{vis}}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}+S_{vis}:\frac{\partial E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\right)d{V}_0}$$ (88)

or

$$H^{elas}={\displaystyle {\int }_{V_0}^{}\left(\frac{\partial E}{\partial \overrightarrow{Y}}:ℭ:\frac{\partial E}{\partial \overrightarrow{Y}}+S_{elas}:\frac{{\partial }^{2}E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\right)d{V}_0}$$ (89)

$$H^{vis}={\displaystyle {\int }_{V_0}^{}\left(\frac{\partial E}{\partial \overrightarrow{Y}}:\mathfrak{N}:\frac{\partial E}{\partial \overrightarrow{Y}}+S_{vis}:\frac{{\partial }^{2}E}{\partial \overrightarrow{Y}\otimes \partial \overrightarrow{Y}}\right)d{V}_0}$$ (90)

in which $S_{elas}$ and $ℭ$ are given by Eqs. (42) and (50), and $S_{vis}$ and $\mathfrak{N}$ are given by Eqs. (58) and (59).

As mentioned after equation (48), the difference between solid and fluids is defined by the assumed elastic and viscous constants, i.e.: with the proposed formulation one is able to model (in an isothermal way) both elastic solids, viscoelastic Kelvin-like solids, and compressive viscous fluids.

4 NUMERICAL EXAMPLES

In this section, we present some selected examples to verify the proposed formulation and to explore some of its possibilities. As fluid results are more easily compared with literature 2D solutions, only one 3D qualitative fluid example is presented. After validating the finite strain viscous behavior in 2D fluid analysis, for solid applications, we use three simple 3D examples.

In most cases we adopt ${\gamma }_1={\gamma }_2=1/2$ for shear viscosity. , In order to measure ${\gamma }_i$ qualitative influence, other parameters are also tested in the 3D fluid example. In all fluid examples $TOL=1\text{x}{10}^{-7}$ is adopted together with a $\Delta t$ sufficiently small to limit the number of iterations to 3. As expected, in fluid applications, limiting the number of iterations proved to produce faster processing time than adopting large time steps. For solid examples the same tolerance is adopted, but the number of iterations is not controlled.

4.1 Dam rupture - fluid

This problem is considered a benchmark to test free surface fluid flow solutions and is used to validate the numerical and physical results of the proposed formulation. It is based on the experimental work of Martin and Motce (1958), reproduced numerically by Nithiarasu (2005) using an ALE fluid formulation. The analyzed problem is a dam initially with width $W$ and height $H$, filled with fluid initially at rest. The dam suffers a subtle disrupt at the right wall (Gate), see Fig. 2. The geometric and physical dimensionles properties are taken as presented by Nithiarasu (2005): $W=0.35$, $H=0.70$, $g=1$, $\rho =1$, $\overline{G}=\mu ={10}^{-3}$. As Nithiarasu (2005) treats the fluid as incompressible we adopted a high value for the bulk modulus ($K=2.15\text{x}{10}^9$) to check the formulation overall behavior.

Fig. 2
- Scheme of dam-rupture of Martin and Motce test (1958): (a) initial position; and (b) flow under open gate

Fig. 3 shows the adopted meshes, with number of nodes and element order. The elements are 3D prismatic with unitary thickness and linear approximation in this direction. Both, vertical and horizontal walls are slip walls and the adopted time step is $2.5\text{x}{10}^{-4}$, no dimensions are given by reference. The analysis is carried out using 6700 time steps with the maximum of 3 iterations at each step.

Fig. 3
- FE meshes for dam-rupture model: (a) 38 cubic elements; (b) 324 cubic elements; (c) 392 cubic elements; and (d) 3004 linear elements.

In a first analysis stage, with the intact reservoir, the water is allowed to conform to meet the initial hydrostatic stress distribution. In the second stage, the right side wall (gate) is instantly removed and the fluid is free to flow. We compare the obtained results to the experimental values of Martin and Motce (1958), since the numerical results of Nithiarasu (2005) are practically coincident with ours, see Fig. 4. The dimensionless time used by the reference to make graphics is obtained by the conversion $t*=t\sqrt{2g/W}$.

Fig. 4
– Relative enlargement of the fluid base along time.

As one can observe, the unstructured cubic approximation meshes present excellent results and practically coincide with the results presented by Nithiarasu (2005). The structured cubic approximation mesh also shows very good results, with a small reduction in total displacement noticed only for the poorest mesh. The unstructured mesh with linear approximation presents an excessive displacement due to element’s distortions. These distortions occur to ensure the nearly incompressible behavior, see Fig. 5. No regularization strategies are used for any of the presented results.

Fig. 5
– Detail of FE mesh for dam-rupture using 3004 linear elements: (a) Initial position, (b) Final position

Fig. 6, presents some snapshots for discretization 324-Cubic, with colors representing displacements, from 0 (red) to the maximum displacement at present instant (blue).

Fig. 6
– Displacement fields of dam-rupture model using 392 cubic elements at different times: (a) $t^{*}=0.896$, (b) $t^{*}=2.391$, (c) $t^{*}=4.183$.

In Nithiarasu (2006), an erratum of Nithiarasu (2005), the author presents a very smooth pressure profile along time that can be reproduced by our formulation if a smaller bulk modulus is employed. Although this study is interested in proving the formulation capable of specific applications, the following reasoning is presented for general interest.

For high bulk modulus our formulation exhibits, at the first time step $(\Delta t=2.5\text{x}{10}^{-4})$, after the instantaneous disrupt, a result (presented in Fig. 7 for discretization 392 – Cubic) that indicates the need of a smaller time step in order to verify the pressure wave propagation after a localized instantaneous pressure reduction. This pressure wave is volumetric and is also called sound wave. In Fig. 8, using a time step of $\Delta t=7.0\text{x}{10}^{-7}$, the pressure wave is detected by exhibiting the first 20 time steps of this analysis, being the first static. Remembering that the pressure wave velocity is calculated as $C_p=\sqrt{K/\rho }=4.63\text{x}{10}^4$ it is necessary approximately 10 time steps to the wave front goes through the horizontal distance. As expected, analyzing Fig. 8 the wave arrives at the left wall about the 11th snapshot, after which the reflections start. One may note that in Fig. 8 the color scale is different for each frame, that is, the pressure level is not the same for each presented time.

Fig. 7
- Stress fields using 392 cubic elements (a) Static result – (b) Dynamic result at first time step ($\Delta t=2.5\text{x}{10}^{-4}$)

Fig. 8
- Pressure wave propagation of dam-rupture model using 392 cubic elements for the first 20 time steps$\left(\Delta t=7.0\text{x}{10}^{-7}\right)$.

Using this small time step and this high bulk modulus the numerical result of pressure wave disperses for the total time necessary to represent the “slow” dam flow, so we relaxed the bulk modulus value to $K=215$ that keeps the low compressibility of the fluid and ensure the stabilized pressure calculation. This procedure is justified by the artificial wave speed (directly related to an artificial compressibility) adopted in Nithiarasu (2005) to capture quasi static pressures in incompressible fluids. In Fig. 9 the achieved pressure profile is qualitatively compared with the one presented by Nithiarasu (2006) as the reference does not presents time or values. The signal of the pressure is negative because it is calculated here as in solid mechanics, i.e., the hydrostatic part of the stress tensor.

Fig. 9
- Comparison of pressure profile of dam-rupture models: (a) Nithiarasu (2006); and (b) proposed model

It should be noted that, since there is no shear wave in the fluid, when the bulk modulus is high, the time intervals required for the proper modeling of pressure waves and the overall movement of viscous fluids are very different. It is also noteworthy that the discretization with linear approximation did not present precise results in our formulation, see the end of example 4.3 for an attempt to a geometrical explanation.

4.2 Solitary wave - fluid

This example is used to validate our formulation in situations of difficult moving boundary situations. In this example, the propagation of a solitary wave in a water tank is studied. The adopted physical properties are (Nithiarasu (2005) and Sung et al. (2000)): $\rho =1000\text{kg}/{\text{m}}^3$, $\overline{G}=\mu =0.001\text{Pa}\text{s}$ and $K=2.15\text{GPa}$. The problem consists of water confined between three smooth walls, with free surface in the top, as shown in Fig. 10. The adopted dimensions for the numerical analysis are (Nithiarasu (2005) and Sung et al. (2000)): $\mathcal{l}=16.0\text{m}$, $H=0.3\text{m}$ and $d=1.0\text{m}$.

Fig. 10
- Scheme of solitary wave model of Laitone (1960)

The tank thickness is $1m$ and the adopted gravity acceleration is $g=9.8\text{m}/{\text{s}}^2$. Initial conditions are calculated based on the analytical solutions given by Laitone (1960) for $t=0\text{s}$, i.e.:

$${\dot{Y}}_1=v_x=\sqrt{gd}\frac{H}{d}{\left(\mathrm{sech}\left(\sqrt{\frac{3H}{4d^3}}x\right)\right)}^2$$ and $${\dot{Y}}_2=v_y=\sqrt{3gd}{\left(\frac{H}{d}\right)}^{3/2}\left(\frac{y}{d}\right){\left(\mathrm{sech}\left(\sqrt{\frac{3H}{4d^3}}x\right)\right)}^2\tanh \left(\sqrt{\frac{3H}{4d^3}}x\right)$$

And introduced directly in equations (78) for $s=0$.

It is important to mention that the initial shape of the solitary wave is given by:

$$Y_2=h=d+H\mathrm{sech}\left(\sqrt{\frac{3H}{4d^3}}(x-ct)\right)$$ with $$c=1+\frac{1}{2}\frac{H}{d}-\frac{3}{20}{\left(\frac{H}{d}\right)}^2$$

In this work, it is prescribed a volumetric force of $b_y=\rho g$ that is responsible to develop the gravity pressure. We adopt $\Delta t=0.01s$ and two FEM meshes with cubic approximation, see Figs. 11a and 11b.

Fig. 11
- FE meshes of solitary wave model: (a) coarse (64 cubic elements), (b) fine (138 cubic elements); and linear (414 linear elements).

Fig. 12 shows the vertical displacement of the upper right and left points for the two discretizations, without considering surface tension. Fig. 13 shows the effect of the water surface tension ($t_s=\mathrm{0,072}N/m$) using the richer mesh.

Fig. 12
- Vertical displacements given by our solution (coincident with Nithiarasu (2005) )

Fig. 13
- Surface tension (St) influence - richer mesh

The solution of the richer discretization, without considering the surface tension, has such a good agreement with the solution presented by Nithiarasu (2005) and Sung et al. (2000) that it was not necessary to provide the solution of the references in the figures. In addition, it is observed that the solution with poorer discretization is already sufficiently precise for this problem.

The application of surface tension is performed in a specific way for two-dimensional simulations, multiplying the surface tension by the semi-thickness of the problem ($1/2\text{m}$ in this case) and applying it as forces on the boundary nodes in the direction of its neighboring node. As expected, due to the dimensions of the problem, the surface tension has little influence, even so, there is a reduction in the peak of the second wave (on the right point) and a slight smoothing of the overall response.

The high-order structured mesh is adopted, because, for the first example, it presents very good results. To stress that linear elements are not recommended to be used with our formulation, a structured mesh with linear approximation is also tested for the solitary wave problem, see Fig 11c.

Fig. 14 shows position snapshots for different instants and different discretizations, including the inadequate linear approximation.

Fig. 14
- Comparison of displacement fields of solitary wave model using different meshes at three times: (a) t=2.4s; (b) t=7s; and (c) t=11.6s.

The presented maximum displacements (at Fig. 14) correspond to the finer cubic discretization are $u_y(2.4\text{s})=0.722\text{m}$, $u_y(7.0\text{s})=0.671\text{m}$, $u_y(11.6s)=0.675m$.

Comments about the bad behavior of linear approximation are given at the end of the next example.

4.3 Formation of a flat bubble – fluid

This example is used to confirm and explain the bad behavior of the linear element regarding our fluid applications (locking). It consists of a $1.0{\text{m}}^{\text{3}}$water volume, initially in the shape of a unit-side cube, released to move in the directions $x_1$ and $x_2$, under the sole action of the considered surface tension $t_s=72\cdot {10}^{-3}\text{N/m}$. The application of this surface tension is made as explained in the previous example and its corresponding nodal force value is $f=36\cdot {10}^{-3}\text{N}$.

The water properties are: density $\rho =1000{\text{kg/m}}^{\text{3}}$, shear viscosity $\overline{G}=\mu =1.02\text{x}{10}^{-3}\text{Pas}$, bulk modulus $K=2.15\text{GPa}$ and volume viscosity $\overline{K}=2.5\text{x}{10}^{-3}\text{Pa}\text{s}$ , see Holmes et al. (2011), for instance. The adopted time interval is $10\text{ms}$.

Two types of analyses are preformed, without any additional damping (for which the bubble vibrates indefinitely) and other introducing an artificial mass proportional damping ($c=1\text{kg}{\text{.(ms)/m}}^{\text{3}}$) in order to find a final static position. In addition, a relaxation in the water compressibility is adopted in order to verify its influence in the poorest cubic discretization, in the linear discretization in the pressure behavior (eliminating pressure wave effects). In the damped analysis, the water volume converges to the circular shape, as shown in Figs. 15 and 16. Two cubic approximations at the plane $(x_1,x_2)$with linear approximation along thickness $x_3$ are adopted. One with 4x4 elements totalizing 338 nodes and other with 8x8 elements and 1250 nodes.

Fig. 15
– Bubble behavior along time $K=2.15\text{GPa}$ (cubic approximation)

Fig. 16
– Position behavior of the analyzed bubble along time using a mesh of 4x4 cubic elements and bulk modulus $K=2.15\text{x}{10}^{-4}\text{GPa}$

As one can see in Fig. 15, when the fluid is almost incompressible, a more refined cubic approximation is necessary to reproduce the expected behavior. Using the 4x4 cubic discretization, the same problem is solved by increasing the compressibility of the fluid, i.e., reducing the bulk modulus to $K=2.15\text{x}{10}^{-4}\text{GPa}$, see Fig. 16.

The pressure for the 8x8 cubic approximation and near incompressible bulk modulus ($K=2.15\text{GPa}$) and for the 4x4 cubic approximation with bulk modulus $K=2.15\text{x}{10}^{-4}\text{GPa}$ oscillates after the displacement convergence, see Fig. 17. This fact occurs because the wave speed of volumetric stresses (sound wave) is high and the time for their mitigation (even using the volumetric viscosity) is also high, not reached by our analysis.

Fig. 17
– Pressure values ($GPa$) at time $t=120\text{ms}$using a mesh of 8x8 cubic element and considering a damping of $c=1\text{kg}{\text{.(ms)/m}}^{\text{3}}$.

Only in the case with large compressibility ($K=2.15\text{x}{10}^{-9}\text{GPa}$), for both, cubic (Fig. 17) and linear (Fig. 18) approximations, the stresses distribution is very close to the expected static value, calculated analytically as $p=-1.44\text{x}{10}^{-10}\text{GPa}$. The signal of the pressure is negative because it is calculated here as in solid mechanics, i.e., the hydrostatic part of the stress tensor.

Fig. 18
– Analysis of the volumetric locking considering various bulk modulus and a mesh of 24x24 linear elements

For the compressible case$K=2.5\text{x}{10}^{-9}\text{GPa}$, the pressure wave has a very low speed and the corresponding behavior is almost static.

Fig. 18 shows the final artificially damped result obtained for a structured linear mesh with 24x24 elements and 1250 nodes. As one can see, the linear element locks completely, even for low fluid compressibility $K=2.15\text{x}{10}^{-6}\text{GPa}$. The behavior mimics a molded truss in the hydrostatic stress distribution, making the continuum rigid. Only when imposing a very compressive behavior ($K=2.15\text{x}{10}^{-9}\text{GPa}=2.15\text{Pa}$) that elements unlock and the results become coherent.

It should be noted that volumetric locking is a known phenomenon that is easily established by a constraint equation that represents the null volumetric variation $\delta J=0$, which for finite element of any order is written simply as:

$$\frac{\partial J}{\partial E}:\delta E=\frac{\partial J}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}=0$$ (a)

as $\delta \overrightarrow{Y}$ is arbitrary, results the following nonlinear system of equations written in the finite element degrees of freedom:

$$\frac{\partial J}{\partial E}:\frac{\partial E}{\partial \overrightarrow{Y}}(\overrightarrow{Y})=\overrightarrow{0}$$ (b)

This system of equations indicates (due to the imposed volumetric restriction) that a degree of freedom is removed from the set of degrees of freedom of the finite element. As an attempt to do a geometrical interpretation to this phenomenon, the identification of this restriction can be done by leaving only one node of the finite element free and studying its movements. In two-dimensional space and high-order elements, this movement is a curve and allows the accommodation of nodes for more or less complex problems, as discretization is increased or compressibility is reduced.

For linear elements, this accommodation is more limited, even when compressibility relaxation is imposed, see Fig. 18. For incompressible or low compressible situations the volumetric locking is easily identified by the 2D representation of Fig 19. Each triangle (with different colors) represents the same black initial linear finite element and, as its volume does not change ($\delta J=0$), the detached node is enforced to move over a straight line. Thus, a node that should have two degrees of freedom has only one due to the linear approximation of finite elements. The same behavior is easily represented for 3D situations if the finite element is a linear tetrahedral.

Fig. 19
– Using a free node constrained movement to illustrate the volumetric locking in a 2D representation

Finally, this numerical example demonstrates that, even for a large bulk modulus reduction, linear elements lock. Note that we use only one integration point in the plane of the linear element and, thus, there is no way to apply the reduced integration technique. From these reasons, in the next examples, the use of linear elements is abandoned.

4.4 Oil passing through a funnel - fluid:

After validating the proposed formulation for viscous compressible flow and describing the choice of high order elements, we solve this free qualitative example in order to show future possibilities of the proposed model. A high viscosity oil $\overline{G}=\mu =0.2356\text{Pa}\text{s}$ with density $\rho =846.6{\text{kg/m}}^{\text{3}}$ and high compressibility $K=2.15\text{x}{10}^4\text{Pa}$ is subject to gravity acceleration $g=10{\text{m/s}}^{\text{2}}$ and is suddenly released to flow through the funnel shown in Fig. 20a. The oil touches the horizontal surface and flows without friction until it stops and starts to return (due to the surface tension). The adopted surface tension is $t_s=0.032\text{N/m}$.

Fig. 20
– Geometry and discretization used to model the oil passing through a funnel, (a) Dimensions and flow direction; and (b) Initial adopted mesh - Cubic - Structured 12x12

Surfaces are frictionless without the possibility of detachment, that is, dry to shear and wet to detachment. No previous static analysis was performed, that is, the gravity occurs suddenly at the beginning of the analysis. To model the fluid only elements with cubic approximation at faces and structured mesh are adopted (see Fig. 20b), since in the first two examples this strategy proved to be sufficiently precise. The vertical and inclined surfaces are modeled by direct constraints and the horizontal surface is modeled using the penalization strategy.

We adopted $\Delta t=0.25\text{ms}$ until the fluid reaches the horizontal surface, i.e., $67.5\text{ms}$, from which $\Delta t=0.125\text{ms}$ is used to increase stability and keep the number of iteration equal to 3, as mentioned in the introduction of this section.

The evolution of the fluid movement is presented in Fig. 21. Pressure values are presented stressing that this fluid has a very high viscosity and not so high bulk modulus, therefore pressure results are stable without using any relaxation. The signal of the pressure is negative because it is calculated here as in solid mechanics, i.e., the hydrostatic part of the stress tensor.

Fig. 21
– Pressure fields at various time (${\text{g/[(ms)}}^{\text{2}}\text{cm]}$) for the oil passing through a funnel model

Fig. 21 shows a smooth response with a large inertial spread of the fluid. In the return response, due to the surface tension, a configuration is reached that indicates a tendency of almost rupture of the surface, which would promote the separation of a drop of fluid at the end of the flow. To model this rupture is beyond the objectives of this work.

4.5 Simplified Slump test - fluid

This example shows the applicability of the proposed formulation in a simple 3D fluid problem, as well as, the model's possibilities in representing fluids whose viscosity depends on the evolution of shear strains. It is a simplified simulation of a slump test of a cement paste whose physical properties are given by Bouvet et al. (2010): $\overline{G}=\mu =1.4\text{Pa}\text{s}$, $\rho =2500{\text{kg/m}}^{\text{3}}$ and $g=10{\text{m/s}}^{\text{2}}$. We adopted a bulk modulus of $K=215\text{kPa}$ that is sufficient to consider near incompressibility. The employed time step is $\Delta t=0.1\text{ms}$ with total analysis time of $t=0.6s$. The test consists of retaining the cement paste in an inverted hollow bucket, with its subsequent instantaneous release on a smooth surface. Fig. 22 depicts the initial configuration geometry ($D_{\inf }=80\text{mm}$, $D_{\sup }=70\text{mm}$ and $H=40\text{mm}$) and the applied discretizations. Due to the symmetry around the vertical axis, only 1/4 of the geometry is discretized by curved elements with cubic approximation in both basis and height. In discretization (a) the base and the top are divided into 8 curved triangular elements with 4 divisions in height, totalizing 32 prismatic finite elements of cubic approximation (637 nodes). In discretization (b) the base and the top are divided into 16 elements with 4 height divisions, totalizing 64 prismatic elements of cubic approximation (1183 nodes). Observing Fig. 23, for $\gamma =0.5$, the maximum difference in diameter at the end of the analysis is less than 1.5% and the maximum difference in height is less than 0.06%. Thus we concluded that using discretization (b) is sufficient to solve the other cases of example.

Fig. 22
– Simplified slump test geometry and discretizations

Fig. 23
– Vertical displacement fields for different meshes and time, $D_{\inf }$ is the diameter of the base

In this example three different exponents for the viscous law are adopted: ${\gamma }_1={\gamma }_2=0.5$ (quasi Newtonian), ${\gamma }_i=0.25$ and ${\gamma }_i=0.75$ (non-Newtonian). Fig. 23 show some top views for selected instants for both discretizations and all viscous parameters. In Fig 23 we indicate the diameter of the base for each selected time step. The color in maps shows the vertical displacement.

The experimental final diameter (at $t=0.2s$) informed by Bouvet et al. (2010) is 169mm, that includes a residual shear stress of ($20Pa$). The diameter achieved by the proposed formulation (${\gamma }_i=0.75$) at the same time $(t=0.2s)$ is near to the experimental result. However, the proposed formulation solution continues to flow as it does not consider residual stresses. For further developments, we think being possible to include residual stresses in the proposed formulation, enabling the appropriate modeling of pastes and gels.

Figure 24 shows the time history of the relation $D_{\inf }^{num}(t)/D_{final}^{exp}$ for the three adopted rheological models, i.e., ${\gamma }_i=0.25$, ${\gamma }_i=0.5$ and ${\gamma }_i=0.75$.

Fig. 24
- Lower diameter behavior along time.

As one can see, the viscous behaviors are quite different, remembering that the case (${\gamma }_i=0.5$) is the quasi Newtonian one.

4.6 Creep test - solid

This example is adopted to show the unified characteristic of the formulation, i.e., solving a solid problem instead of a fluid problem using the same approach. In this example, a 3D block of unitary dimension is subjected, at its free face, to a constant tensile force in $z=x_3$ direction, see Fig. 25 $\left(F=40\text{kN}\right)$ and its longitudinal strain (${\lambda }_3-1$) is evaluated along time (creep test). Figure 25 presents the test geometry and the discretization of only $1/4$ of the specimen to take advantage of symmetry. We used 4 prismatic elements with linear approximation totalizing 12 nodes. The following physical parameters are adopted for all analyses: $K=1.5\text{MPa}$, $G=9\text{KPa}$, $\rho =0$, $g=0$, $\overline{K}=0$ and ${\gamma }_1={\gamma }_2=1/2$. We used 3 different values for shear viscosity: ${\overline{G}}_a=2\text{KPas}$, ${\overline{G}}_b=1.1\text{KPas}$ and ${\overline{G}}_c=0.2\text{KPas}$ and the corresponding time intervals: $\Delta t_a=0.4\text{s}$, $\Delta t_b=0.2\text{s}$ and $\Delta t_c=0.04\text{s}$. The adopted tolerance is $Tol={10}^{-7}$.

Fig. 25
– Creep tesr schematic representation and adopted discretization (linear elements).

Figure 26 presents the time response of the stretch in z direction and snapshots of the Cauchy viscous stresses ${\sigma }_{11}^{vis}$ and ${\sigma }_{33}^{vis}$ for case (c). As expected, one can see that both viscous stresses ${\sigma }_{11}^{vis}$ and ${\sigma }_{33}^{vis}$ approaches zero at the end of the analysis.

Fig. 26
- Time response of stretch and viscous stress snapshots.

In Fig. 27 one can see the time response of viscous stress ${\sigma }_{33}^{vis}$ and elastic stress ${\sigma }_{33}^{elas}$, both in Cauchy space. It is worth noting that the final cross section continues to be a square with sides $0.360\text{m}$ (full dimension) and corresponds to the final achieved Cauchy stress value. As one can see the proposed formulation is capable of modeling the viscoelastic behavior of highly deformable bodies that do not present an instantaneous response (Kelvin-like model). It is interesting to mention that the strain levels reached by the material are quite high.

Fig. 27
- Viscous and elastic Cauchy stresses along time (${\sigma }_{33}$).

4.7 Indirect relaxation test - solid

In this example, the specimen of the previous example is made up of two materials, see Fig. 28. The basic boundary conditions can be seen in Fig. 25a. Material 1 is viscoelastic with the following properties: $K=15\text{MPa}$, $G=9\text{KPa}$, ${\gamma }_1={\gamma }_2=1/2$ and $\overline{G}=200\text{KPas}$. Material 2 is elastic with the following properties $K=0.15\text{MPa}$ and $G=9\text{KPa}$. In this analysis we adopted $\rho =0$, $g=0$ and $\overline{K}=0$. Figure 28 also illustrates the adopted discretization, i.e., 4 prismatic elements with cubic approximation (112 nodes).

Fig. 28
– Indirect relaxation test: (a) Material distribution; and (b) Discretization for a quarter of the problem (Symmetry)

Initially, material 1 is maintained with all nodes restricted and the free face of material 2 (elastic material) is stretched from coordinate $z=1\text{m}$ through coordinate $z=3.4\text{m}$ using 120 equally spaced steps. Figure 29 shows the achieved configuration with the Cauchy elastic stress distribution ${\sigma }_{33}^{ela}$. Note that as all nodes of material 1 are constrained in the initial position and no stresses are present, but an average value appears at the interface nodes in Fig. 29a due to post processing. Figure 29b presents the viscous and elastic Cauchy stresses at point $(0.167;0.167;0.167)\text{m}$ of material 1 and the elastic Cauchy stress at point $(0.167;0.167;0.833)\text{m}$ of material 2. These coordinates correspond to the initial configuration.

Fig. 29
- Indirect relaxation test: (a) Initial position; and (b) Cauchy stress along time

As one can see in Fig. 29b, the stress developed in material 2 represents the relaxation behavior of the elastic part of the studied specimen. After the total relaxation, the elastic stresses (or total stress, as there is no viscous stress at the end of the analysis) in both materials approximate each other. It is important to stress that the cross section of material 1 varies from $1{\text{m}}^{\text{2}}$ through $0.31{\text{m}}^{\text{2}}$ which explains the large difference between the initial viscous stress and the final elastic stress of material 1. The adopted tolerance is $Tol={10}^{-7}$ and the time interval is $\Delta t=0.4\text{s}$.

4.8 Viscoelastic sandwich circular plate - solid:

In this example we show the use of viscoelastic properties of a core to describe the induced damping in structural elements. Differently from usual structural analysis the damping is stablyshed directly from reological material properties. A simple supported circular plate with radius $R=1\text{m}$ and thickness $t=3\text{cm}$ is subjected to a transverse uniform loading. Only $1/4$ of the structure is modeled using 300 prismatic finite elements with cubic approximation parallel to the plate surface and linear along thickness, totalizing 3 unitary layers, see Fig. 30. The simple support condition is applied at boundary nodes of the bottom face. The load is applied as a volume force on the superior layer of the plate. When viscosity is considered we used ${\gamma }_1={\gamma }_2=0.5$, i.e., Kelvin-Voigt-like viscoelastic model. Three situations are considered:

Fig. 30
– Elastic case: Boundary conditions and transverse displacement $(w_{\max }=0.7044\text{cm})$

1. i

   Only to verify the discretization, the three layers are considered elastic (steel) with properties: $E=200\text{GPa}$ and $\nu =0.25$ that corresponds to $K=133\text{GPa}$ and $G=80GPa$. The adopted transverse load is $b_3=5000{\text{kN/m}}^{\text{3}}$ on the top lamina, corresponding to a surface force of $h_3=50{\text{kN/m}}^{\text{2}}$. For this case the achieved central transverse displacement is $w=0.7043\text{cm}$, $2.9\%$ larger than the Kirchhoff kinematics analytical solution that is $w_k=0.684\text{cm}$ see Timoshenko and Woinowsky-Krieger, S. (1959) for instance.

This result is expected as the adopted solid element is more flexible than the Kirchhoff kinematics. It is important to mention that this element has ben validated for solid analysis in Carrazedo et al. (2020) and Carrazedo and Coda (2018).

1. ii

   Keeping the loading of case (i), the material of the central layer is substituted by Polypropylene with the elastic properties $E=1.088\text{GPa}$ and $\nu =0.49$ that correspond to $K=18.13\text{GPa}$ and $G=0.365\text{GPa}$. The adopted shear viscosity property is $\overline{G}=6.756\text{GPa×s}$. These values are adapted from Fazekas and Goda (2018). The central displacement for the elastic and viscoelastic cases are shown in Fig 31, being the maximum elastic displacement $w_{\max }=0.7432\text{cm}$, $5.52\%$ larger then the steel case (i). We used 100 time steps of $\Delta t=0.01\text{s}$ for the viscoelastic analysis.

   

   Fig. 31
   – Sandwich viscoelastic plate: (a) Final displacement field, (b) Viscoelastic displacement at the plate center.
   

1. iii

   For the material of case (ii), considering the steel density ${\rho }_{steel}=7000{\text{kg/m}}^{\text{3}}$ and the Polypropylene density ${\rho }_{pol}=910{\text{kg/m}}^{\text{3}}$, we perform a dynamic analysis considering the same load of previous cases suddenly applied. The central displacement along time is compared with the elastic result of case (ii) at Fig. 32. We used 500 time steps of $\Delta t=0.01\text{s}$.

   

   Fig. 32
   – Central displacement along time for the sandwich plate viscoelastodynamic analysis
   

As one can see by this example, the formulation is capable to simulate viscoelastic (Kelvin-Voigt-like model) structural elements including dynamics.

6 CONCLUSIONS

This work presents a unified constitutive model for both low compressible isothermal viscous fluids and simple Kelvin-Voigt-like viscoelastic solids. Interpreting the Lagrangian strain directions - derived from the Flory’s decomposition for hyperelastic materials - we proposed a coherent dissipative virtual work for the internal energy allowing viscous analysis. Combining the hyperelastic potential and the dissipative internal energy, both elastic and viscous stresses are calculated in a straightforward way. The proposed constitutive model is successfully implemented in the positional finite element method, resulting in a total Lagrangian formulation whose nodal parameters are positions. Two-dimensional benchmarks from literature are used to verify the proposed constitutive model and the developed numerical formulation to represent viscous fluids. A flat bubble example is used to show the locking relief of the proposed numerical formulation by adopting high order approximations. A simple slump test is used to show the applicability of the model to 3D compatible problems. Static and dynamic viscoelastic 3D solid examples are used to demonstrate the unified characteristic of the proposed approach. As the formulation is based on scalar quantities (Lagrangian strain invariants), the proposed constitutive model opens the possibility of furthers developments to include residual stresses to model pastes, clay and gels.

ACKNOWLEDGEMENTS

This study has a major support grant #2020/05393-4 of the São Paulo Research Foundation (FAPESP) and has a graduate grant from the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

References

• Akkerman, I., Bazilevs, Y., Calo. V.M., et al (2008). The role of continuity in residual-based variational multiscale modeling of turbulence. Computational Mechanics, v. 41 (3), p. 371-378.
• Anderson, J. D. (1995). Computational fluid dynamic - the basics with applications, McGraw-Hil, New York.
• Argyris, J.H. (1954). Energy theorems and structural analysis Part 1. Aircraft Eng., v. 26, p. 383
• Argyris, J., Doltsinis, I.St., Silva V.D. (1991). Constitutive modelling and computation of non linear viscoelastic solids. Part I: Rheological models and integration techniques, Comput. Methods Appl. Mech. Engrg. v. 88, p, 135–163.
• Arruda, E.M.; Boyce, M.C. (1993). A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids, v. 41 (2), p. 389-412.
• Benson, D.J., Bazilevs, Y., Hsu, M.C., et al. (2011). A large deformation, rotation-free, isogeometric shell, Computer Methods in Applied Mechanics and Engineering, v. 200 (13-16), p. 1367-1378
• Bischoff, M., Ramm, E. (2000). On the physical significance of higher order kinematic and static variables in a three-dimensional shell formulation, International J. Solids Struct., v. 37, p. 6933-6960
• Bonet, J., Wood, R. (1997). Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge University Press, New York.
• Bouvet, A., Ghorbel, E., Bennacer, R. (2010). The mini-conical slump flow test: Analysis and numerical study, Cement and Concrete Research, v. 40, p. 1517 -1523.
• Brooks, A.N., Hughes, T.J. (1982). Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering, v.32, p. 199–259.
• Carrazedo, R., Coda, H.B. (2018). Triangular based prismatic finite element for the analysis of orthotropic laminated beams, plates and shells, Composite Structures‏, v. 168. P. 234-246
• Carrazedo, R., Paccola, R.R., Coda, H.B. (2020). Vibration and stress analysis of orthotropic laminated panels by active face prismatic finite element, Composite Structures‏, v. 244, e. 112254
• Chen, W.H., Chang, C.M., Yeh, J.T. (1993). An incremental relaxation finite element analysis of viscoelastic problems with contact and friction, Comput. Methods Appl. Mech. Engrg., v. 9, p. 315–319
• Chung, W.J., Cho, J.W., Belytschko, T. (1998). On the dynamic effects of explicit FEM in sheet metal forming analysis, Engineering Computations, v. 15 (6-7), p. 750-776
• Chung, T.J. (2002). Computional fluid dynamics, Cambridge University Press, Cambridge, UK.
• Clough, R.W. (1960). The Finite Element Method, in Plane Stress Analysis, Proc. 2nd A.S.C.E. Conf: on Electronic Comp., Pittsburgh, PA.
• Clough, R.W. (1989). Original formulation of the Finite Element Method, Proc. ASCE Structures Congress Session on Computer Utilization in Structural Eng. San Francisco, pp 1-10.
• Coda, H. B. (2003). An exact FEM geometric non-linear analysis of frames based on position description. In: 17th International conference of mechanical engineering, São Paulo, Brazil.
• Coda, H.B., Greco, M. (2006). Positional FEM formulation for flexible multi-body dynamic analysis, Journal of Sound and Vibration, v. 290 (3-5), p. 1141-1174
• Coda, H.B., Paccola, R.R. (2008). A positional FEM Formulation for geometrical non-linear analysis of shells, Latin American Journal of Solids and Structures, v. 5 (3), p. 205-223
• Coda, H.B. (2009). Two dimensional analysis of inflatable structures by the positional FEM, Latin American Journal of Solids and Structures, v. 6 (3), p. 187-212
• Coda, H.B., Paccola, R.R. (2011), A FEM procedure based on positions and unconstrained vectors applied to non-linear dynamic of 3D frames‎, Finite Elements in Analysis and Design, v. 47 (4), p.‏ 319-333.
• Coda, H.B. (2015). Continuous inter-laminar stresses for regular and inverse geometrically non linear dynamic and static analyses of laminated plates and shells., Composite Structures, v. 132, p. 406-422.
• Siqueira, T.M. and Coda, H.B. (2019). Flexible actuator finite element applied to spatial mechanisms by a finite deformation dynamic formulation, Computationa Mechanics, v. 64 (6), p. 1517-1535
• Coda, H.B. (2021). An alternative finite strain elastoplastic model applied to soft core sandwich panels simulation, Latin American Journal of Solids and Structuresv. 18 (6), e392
• Courant, R. (1942). Variational methods for the solution of problems of equilibrium and vibrations, Trans. Amer. Math. Soc. P. 1-23
• Donea, J., Giuliani, S., Halleux, J.P. (1982). An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions, Computer Methods in Applied Mechanics and Engineering, v. 33 (1-3), p. 689–723.
• Duarte. F., Gormaz. R., Srinivasan. N. (2004). Arbitrary Lagrangian-Eulerian method for Navier-Stokes equations with moving boundaries, Comput. Methods Appl. Mech. Eng., v. 193, p. 4819–4836
• Düster, A., Hartmann, S., Rank, E. (2003). p-FEM applied to finite isotropic hyperelastic bodies, Comput. Methods Appl. Mech. Eng., v. 192, p. 5147–5166.
• Elias, R.N., Martins, M.A.D., Coutinho, A.L.G.A. (2006). Parallel edge-based solution of viscoplastic flows with the SUPG/PSPG formulation, Computational Mechanics, v. 38 (4-5), p. 365-381.
• Fazekas, B., Goda, T. (2018). Characterisation of large strain viscoelastic properties of polymers. Bánki Közlemények, v. 1 (1), Hungary.
• Flory, P.J. (1961). Thermodynamic relations for high elastic materials. Transactions of the Faraday Society, v.57, p. 829-838
• Franci, A., Oñate, E., Carbonell, J.M. (2016). Unified Lagrangian formulation for solid and fluid mechanics and fsi problems, Comput. Methods Appl. Mech. Eng., v. 298, p. 520–547.
• Gasser, T.C., Holzapfel, G.A. (2002). A rate-independent elastoplastic constitutive model for biological fiber-reinforced composites at finite strains: continuum basis, algorithmic formulation and finite element implementation, Computational Mechanics v. 29 (4-5), p. 340-360.
• Gent, A.N. (2001).Engineering With Rubber: How to Design Rubber Com-ponents", Hanser Gardner, 2nd ed.
• Gruttmann, F., Wagner, W. (2001). Shear correction factors in Timoshenko's beam theory for arbitrary shaped cross-sections, Computational Mechanics, v. 27 (3), p. 199-207.
• Hartmann, S., Neff, P. (2003). Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility. International Journal of Solids and Structures, v. 40 (11), p. 2767-2791.
• Havner, K.S. (1966). On formulation and iterative solution of small strain plasticity problems, Quarterly of Applied Mathematics, v. 23 (4), p. 323-335.
• Holmes, M.J. et al. (2011). Temperature dependence of bulk viscosity in water using acoustic spectroscopy, J. Phys.: Conf. Ser. V. 269.
• Holzapfel, G.A .(1996). On large strain viscoelasticity: continuum formulation and finite element applications to elastomeric structures, Int. J. Numer. Methods Eng., v.39, p. 3903–3926.
• Holzapfel, G.A., Simo, J.C. (1996). Entropy elasticity of isotropic rubber-like solids at finite strains, Computer Methods in Applied Mechanics and Engineering, v. 132 (1-2), p. 17-44.
• Hudobivnik, B., Aldakheel, F., Wriggersm P. (2019). A low order 3D virtual element formulation for finite elasto-plastic deformations, Computational Mechanics, v. 63 (2), p. 253-269.
• Idelsohn, S.R., Oñate, E., Pin, F.D. (2004). The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves, Int. J. Numer. Meth. Engrg., v. 61, p. 964–989.
• Idelsohn, S.R., Oñate, E., Pin, F.D., Calvo, N. (2006). Fluid-structure interaction using the particle finite element method, Comput. Methods Appl. Mech. Eng., v. 195, p. 2100–2123.
• Idelsohn, S.R., Marti, J., Limache, A., Oñate, E. (2008). Unified Lagrangian formulation for elastic solids and incompressible fluids: Application to fluid–structure interaction problems via the PFEM, Comput. Methods Appl. Mech. Engrg., v. 197, p. 1792–1776.
• Jiao, Y., Fish, J. (2018). On the equivalence between the multiplicative hyper-elasto-plasticity and the additive hypo-elasto-plasticity based on the modified kinetic logarithmic stress rate, Computer Methods in Applied Mechanics and Engineering, v. 340, p. 824-863.
• Kröner, E. (1959). Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Arch. Ration. Mech. Anal. v. 4, p. 273.
• Laitone, E.V. (1960). The second approximation to cnoidal and solitary waves, Journal of fluid mechanics, v. 9, p. 430-444.
• Lánczos, C. (1970). The Variational Principles of Mechanics, Dover, New York USA.
• Latorre, M., Montans, J.F. (2015) Anisotropic finite strain viscoelasticity based on the Sidoroff multiplicative decomposition and logarithmic strains, Computational Mechanics, v. 56 (3), p. 503-531.
• Lee, E.H.(1969). Elastic–plastic deformations at finite strains. Journal of Applied Mechanics (ASME), v. 36, p. 1– 6.
• Martin, J.C., Motce, W.J. (1958). An experimental study of the collapse of liquid columns on a rigid horizontal plane, Philosophical Transactions of the Royal Society of London, Series A 244: 312–324.
• Mesquita, A.D. and Coda, H.B. (2007a). A boundary element methodology for viscoelastic analysis: Part I with cells, Applied Mathematical Modelling v. 31 (6), p.1149-1170
• Mesquita, A.D. and Coda, H.B. (2007b). A boundary element methodology for viscoelastic analysis: Part II without cells, Applied Mathematical Modelling v.31 (6), p.1171-1185
• Miehe, C., Aldakheel, F., Mauthe, S. (2013). Mixed variational principles and robust finite element implementations of gradient plasticity at small strains, International Journal for Numerical Methods in Engineering, v. 94 (11), p. 1037-1074.
• Mooney, M., (1940). A theory of large elastic deformation. Journal of Applied Physics, v. 11 (9), p. 582-592.
• Nithiarasu, P. (2005). An arbitrary Lagrangian eulerian (ale) formulation for free surface flows using the characteristic-based split (cbs) scheme, Int. J. Numer. Meth. Fluids, v. 48, p. 1415–1428.
• Nithiarasu, P. (2006). Erratun An Arbitrary Lagrangian Eulerian (ALE) formulation for free surface ows using the Characteristic Based Split (CBS) scheme (Int. J. Numer. Meth. Fluids 2005; 48:1415–1428), Int. J. Numer. Meth. Fluids, v. 50, p. 1119–1120.
• Ogden, R. W. (1972). Large deformation isotropic elasticity: on the correlation of theory and experiment for compressible rubberlike solids. Proceedings of the Royal Society of London, v. 328, n. 1575, p. 567-583.
• Ogden, R.W. (1984). Non-linear Elastic Deformations. Ellis Horwood, New York.
• Pascon, J.P., Coda, H.B. (2013). High-order tetrahedral finite elements applied to large deformation analysis of functionally graded rubber-like materials,‎ Applied Mathematical Modelling, v. 37 (20-21), p. 8757-8775.
• Pascon, J.P., Coda, H.B. (2015). Large deformation analysis of functionally graded elastoplastic materials via solid tetrahedral finite elements, Computers & Structures, v. 146, p.59-75
• Pascon, J.P., Coda, H.B. (2017), Finite deformation analysis of visco-hyperelastic materials via solid tetrahedral finite elements, Finite Elem. Anal. Des., v. 133, p. 25-41
• Reese, S., Wriggers, P. (1997). A material model for rubber-like polymers exhibiting plastic deformation: Computational aspects and a comparison with experimental results, Computer Methods in Applied Mechanics and Engineering, v. 148 (3-4), p. 279-298.
• Reddy, J.N., Gartling, D.K. (2010). The Finite Element Method in Heat Transfer and Fluid Dynamics. CRC Press, Boca Raton, USA.
• Rivlin, R. S. (1948a). Large elastic deformations of isotropic materials. I. Fundamental concepts. Philosophical Transactions of the Royal Society of London, v. 240, n. 822, p. 459-490.
• Rivlin, R. S. (1948b). Large elastic deformations of isotropic materials. IV. Further developments of the general theory. Philosophical Transactions of the Royal Society of London, v. 241, n. 835, p. 379-397.
• Rivlin, R., Saunders, D. (1951). Large elastic deformations of isotropic materials VII. Experiments on the deformation of rubber, Philos.Trans. R. Soc. London Ser. A, v. 243, p. 251–288
• Sanches, R.A.K., Coda, H.B. (2013). Unconstrained vector nonlinear dynamic shell formulation applied to Fluid Structure Interaction, Computer Methods in Applied Mechanics and Engineering, v. 259, p. 177-196.
• Sansour, C., Bednarczyk, H. (1995). The Cosserat surface as a shell-model, Theory and Finite-Element Formulation, Computer Methods in Applied Mechanics and Engineering, v. 120 (1-2), p. 1-32
• Schwarze, M., Vladimirov, I.N., Reesem S. (2011). Sheet metal forming and springback simulation by means of a new reduced integration solid-shell finite element technology, Computer Methods in Applied Mechanics and Engineering, v. 200 (5-8), p. 454-476.
• Shutov, A.V., Landgraf, R., Ihlemann, J. (2013). An explicit solution for implicit time stepping in multiplicative finite strain viscoelasticity, Computer Methods in Applied Mechanics and Engineering, v. 265, p. 213-225.
• Simo, J.C. (1992). Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, Computer Methods in Applied Mechanics and Engineering, v. 99 (1), p. 61-112.
• Sung, J., Choi, H.G., Yoo, J.Y. (2000). Time-accurate computation of unsteady free surface flows using an ALE-segregated equal-order FEM, Computer methods in applied mechanics and Engineering, v. 190 (11-12), p. 1425-1440.
• Tezduyar, T.E. (1992). Stabilized finite element formulations for incompressible flow computations, Advances in Applied Mechanics, v. 28, p.1-44.
• Tezduyar, T.E., Behr, M., Liou, J. (1992). A new strategy for finite element computations involving moving boundaries and interfaces - the deforming-spatial-domain/space-time procedure: I. The concept and the preliminary numerical tests, Comput. Methods Appl. Mech. Engrg., v. 94, p. 339–351.
• Tezduyar, T.E., Senga, M. (2006). Stabilization and shock-capturing parameters in SUPG formulation of compressible flows, Computer Methods in Applied Mechanics and Engineering, v. 195 (13-16), p. 1621–1632.
• Timoshenko, S.P. and Woinowsky-Krieger, S. (1959). Theory of Plates and Shells, 2nd edn,McGraw-Hill, New York.
• Treloar, L.R.G. (1943). The elasticity of a network of long-chain molecules. I. Transactions of the Faraday Society, v. 39, p. 36-41.
• Turner, M.J., Clough, R.W., Martin, H.C., Topp, L.T. (1956). Stiffness and deflection analysis of complex structures, J. Aeronaut. Scs, v. 25, p. 805-823
• Vergori, L., Destrade, M., McGarry, P., et al. (2013). On anisotropic elasticity and questions concerning its Finite Element implementation, Computational Mechanics, v. 52 (5), p. 1185-1197.
• Yeoh, O.H., (1990). Characterization of elastic properties of carbon-black-filled rubber vulcanizates. Rubber Chemistry and Technology, v. 63, n. 5, p. 792-805.
• Zienkiewicz, O.C., Cheung, Y.K. (1965). The finite elements in the solution of field problems, Engineer, v. 220, p. 507–510.
• Zienkiewicz, O.C, Taylor, R.L., Nithiarasu, P. (2005). The Finite Element Method: Fluid Dynamics, Butterworth Heinemann Linacre house, Oxford, UK.

Edited by

Editor: Rogério José Marczak

Publication Dates

• Publication in this collection
  27 June 2022
• Date of issue
  2022

History

• Received
  22 Feb 2022
• Reviewed
  26 May 2022
• Accepted
  27 May 2022