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

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 β 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.


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 Latin American Journal of Solids and Structures, 2022, 19(4), e449 3/37 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 β 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.

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.

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: Latin American Journal of Solids and Structures, 2022, 19(4), e449 4/37 in which i y is a position of a particle inside the continuum, ij σ is the Cauchy stress component acting on plane i in the direction j , ρ is the mass density, i y  is the particle acceleration following direction i (material derivative of velocity) and i b 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, ( ) Integrating Eq. ( 3) over the spatial domain, the virtual work is written as: and, by distributing the integral, one gets: In order to deal with the last term of Eq. ( 5), we use the following property: Replacing Eq. ( 6) in the last term of Eq. ( 5), and applying the divergence theorem, one writes: where j n is the j th − component of the normal unity vector to the current surface S of the analyzed continuum.
Recalling the Cauchy formula, i ji j h n σ = (in which i h is the distributed force over the Neumann boundary), and taking into account the stress tensor symmetry, we write:  (Ogden(1984) ).Grouping terms of Eq. ( 8) one has: in which K is the kinetic energy, P is the external force potential and Ψ 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.
Latin American Journal of Solids and Structures, 2022, 19(4), e449 5/37 To obtain the Lagrangian description of Equation ( 8), the continuity theorem is applied to the first term of Eq. ( 8), and volume forces are considered proportional to density, In equations ( 10) and ( 11) we use the well known property δ (energetically conjugate to the second Piola-Kirchhoff stress), are related by (Ogden(1984) and Bonet and Wood (1997)): A is the deformation gradient as usually given in finite strain elasticity, J is its determinant and the Green strain is defined as From Equations ( 11) and ( 12), the week form of the equilibrium equation in the Lagrangian description is written as: 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.

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: : 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: ψ and 2 iso ψ are written, respectively, as a function of the first and second invariants of the isochoric part of the Cauchy-Green stretch tensor.The volumetric Ĉ and isochoric 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: with ( ) ( ) 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: results in the desired multiplicative decomposition applied to the Cauchy-Green stretch tensor as: One may observe that 2 det( ) det( ) From C and Ĉ 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: from which one calculates the corresponding second Piola-Kirchhoff stress as: E a symmetric tensor of the same order of the Green strain.It is simple to show that: Applying the first part of Eq. ( 12) over the second Piola-Kirchhoff stress of Eq. ( 23) (pulling forward) results in: ( ) from which one calculates the corresponding second Piola-Kirchhoff stress as: where β is a scalar and 1 iso E is a symmetric tensor of the same order of the Green strain.Opening 1 I one writes: from which one finds: Considering Eq. ( 29) and applying the transformation (12) over Eq. ( 27) one finds: ( ) As ( ) Latin American Journal of Solids and Structures, 2022, 19(4), e449  8/37 ( ) is the first isochoric Lagrangian strain direction corresponding to a deviatory strain, and ( ) 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 2 I as: ( ) from which we calculate the corresponding second Piola-Kirchhoff stress as: with γ and 2 iso E being respectively a scalar and a symmetric tensor of the same order as the Green-Lagrange strain tensor.It is interesting to write: in which 2 I is the second invariant of the Cauchy stretch.With some algebraic effort over Eq. ( 34), results: Using Eq. ( 12) over Eq. ( 33) and considering Eq. ( 35), one writes: From Equation (36), using Eq. ( 37), one calculates, ( ) ( ) Latin American Journal of Solids and Structures, 2022, 19(4), e449  9/37 In equation ( 39) one concludes that 2 iso σ is deviatory and, thus is the second isochoric Lagrangian strain direction, and ( ) is the intensity of the second deviatory Piola-Kirchhoff stress.Thus, 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.: ( ) and In which K is the bulk modulus, G is the transversal elastic modulus and 0 n ≠ is an integer.Other hyperelastic potentials that respect the volumetric and isochoric split can also be chosen.

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.: 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.: Latin American Journal of Solids and Structures, 2022, 19(4), e449 10/37 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, * vis S 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., / : From Eq. ( 46) and Eqs. ( 24), ( 29) and ( 35) the following expression for the second Piola-Kirchhoff viscous stress is written, in which K is the fluid volumetric viscosity Holmes et al. (2011) and i G are shear (isochoric directions 1 and 2) viscosities and vol E , 1 iso E and 2 iso E are strain directions.In order to be coherent with the elasticity understanding we assume . Adopting the viscous parameters 1 2 1/ 2 γ γ = = and 0 K = 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 α and 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 α and i γ is out of this study objectives; however some inspections are made in example 4.5.
. For viscoelastic solids all elastic and viscous constants are taken into account.Moreover, it is important to note that the volumetric viscosity is disregarded 0 K = 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: elas vis elas vis The first term can be calculated as: ∂ ∂E is given by Eq. ( 24), 1 / I ∂ ∂E is given by Eq. ( 29), 2 / I ∂ ∂E is given by Eq. ( 35) and the other terms are given, as: ( ) One may note that in this model ( ) ( ) 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: Differentiating Eq. ( 58) regarding current strain (implicit method), it results the following numeric tangent constitutive tensor at the current time ( 1 s + ): Latin American Journal of Solids and Structures, 2022, 19(4), e449 12/37 Notice that terms related to the past (of Eq. ( 58)) do not appear in Eq. ( 59) as they are constant.Substituting vol E , 1 iso

E and
2 iso E by their respective definitions ( 24), ( 29) and ( 35), and using index notation, it results: in which, for simplicity, the current time is not indicated.

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 Y  as the analysis parameters, thus the continuum position and acceleration are represented by: in which φ are usual FEM shape functions.Volume and surface forces, if desired, can be approximated by their nodal values, as: Latin American Journal of Solids and Structures, 2022, 19(4), e449 13/37 where ϕ are surface shape functions.
The strain variation δ E should also be given as a function of position variation Y δ  , as: Introducing such changes and considereing Eq. ( 48), the weak form of the equilibrium Eq. ( 13) becomes: ( ) As the nodal position variation Y δ  is arbitrary, Eq. ( 66) turns into the following nonlinear system of equations: 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 ext F  .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 ( 1 2 , ξ ξ ), and any order of approximation along its thickness direction ( 3 ξ ).Fig. 1 shows a finite element with linear approximation at 3 ξ direction and cubic approximation at 1 ξ and 2 ξ directions.( ) ( ) in which the last one is an index version of Eq. ( 61).The gradients of these mappings are given by: Thus, the gradient A of the deformation function f  (see Equations ( 12) and ( 16)) is given by: ( ) Eq. ( 70) is numerically written at each integration point for a trial position Y  covering all nodes of the problem.As the formulation is total Lagrangian, one may observe that the gradient of the initial configuration mapping 0 A is calculated only once in the numerical process.Moreover, it is clear the dependence of A regarding Y  .
After calculating the deformation gradient A , one calculates the Cauchy-Green stretch tensor and the Green strain as: ( ) Using the equations of the previous sections, all necessary variables, such as vol E , can be numerically calculated at each integration point of Eq. (66) or Eq. ( 67) and are function of nodal positions 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: in which 1 s t + is the current instant.Thus, one writes Eq. (67) as: ( ) 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: ( ) in which β and γ are free parameters, generally adopted as and 1 / 2 γ = .Isolating the current velocity and the current acceleration from Eqs. ( 74) and ( 75), results: with the following auxiliary values: ( ) Substituting Eqs. ( 76) and ( 77) into ( 73), results: ( ) Eq. ( 79) is understood simply as ( ) , revealing the nonlinear behavior of equilibrium equations regarding ( ) , 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: where ( ) with elas H and din H being, respectively, the elastic and dynamic parts of the Hessian matrix.The viscous term vis H 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 Y ∆  is achieved as: ( ) ( ) At the beginning of each iteration, acceleration and velocity should be recalculated as: are values of the past, they remain unaltered during iterations.The stopping test is given by: in which TOL is the prescribed numerical tolerance.
At the first time step, the initial acceleration is calculated by the dynamic equilibrium equation as: 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: 0 0 2 0 0 : : : : or 0 2 0 : : : in which elas S and C are given by Eqs. ( 42) and ( 50), and vis S and 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.

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 Latin American Journal of Solids and Structures, 2022, 19(4), e449 17/37 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 1 2 1/ 2 γ γ = = for shear viscosity., In order to measure i γ qualitative influence, other parameters are also tested in the 3D fluid example.In all fluid examples 7 1x10 TOL − = is adopted together with a 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.

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)  ) to check the formulation overall behavior.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 4 2.5x10 − , no dimensions are given by reference.The analysis is carried out using 6700 time steps with the maximum of 3 iterations at each step.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 * 2 / t t g W = .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.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 , 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 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.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 215 K = that keeps the low Latin American Journal of Solids and Structures, 2022, 19(4), e449 20/37 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.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.

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)  .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)   The tank thickness is 1m and the adopted gravity acceleration is 2 9.8m / s g = .Initial conditions are calculated based on the analytical solutions given by Laitone (1960) for 0s t = , i.e.: And introduced directly in equations (78) for 0 s = .
It is important to mention that the initial shape of the solitary wave is given by: Latin American Journal of Solids and Structures, 2022, 19(4), e449 21/37 In this work, it is prescribed a volumetric force of y b g ρ = that is responsible to develop the gravity pressure.We adopt 0.01 t s ∆ = and two FEM meshes with cubic approximation, see Figs. 11a and 11b.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.and Structures, 2022, 19(4), e449 22/37 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 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.
The presented maximum displacements (at Fig. 14) correspond to the finer cubic discretization are (2.4s) 0.722m Comments about the bad behavior of linear approximation are given at the end of the next example.

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 3 1.0m water volume, initially in the shape of a unit-side cube, released to move in the directions 1 x and 2 x , under the sole action of the considered surface tension ) 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 1 2 ( , ) x x with linear approximation along thickness 3 x are adopted.One with 4x4 elements totalizing 338 nodes and other with 8x8 elements and 1250 nodes.
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

2.15GPa K =
) and for the 4x4 cubic approximation with bulk modulus 4 2.15x10 GPa K − = 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.
Latin American Journal of Solids and Structures, 2022, 19(4), e449 23/37 ), for both, cubic (Fig. 17) and linear (Fig. 18) approximations, the stresses distribution is very close to the expected static value, calculated analytically as .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.
For the compressible case  ) that elements unlock and the results become coherent.
Latin American Journal of Solids and Structures, 2022, 19(4), e449 25/37 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 0 J δ = , which for finite element of any order is written simply as: as Y δ  is arbitrary, results the following nonlinear system of equations written in the finite element degrees of freedom: 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 ( 0 J δ = ), 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.
Latin American Journal of Solids and Structures, 2022, 19(4), e449 26/37 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.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

0.25ms t ∆ =
until the fluid reaches the horizontal surface, i.e., 67.5ms , from which 0.125ms t ∆ = 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 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.

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 , 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.
Latin American Journal of Solids and Structures, 2022, 19(4), e449 28/37 ) at the same time ( 0.2 ) t s = 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.

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     In this analysis we adopted 0 ρ = , 0 g = and 0 K = .Figure 28 also illustrates the adopted discretization, i.e., 4 prismatic elements with cubic approximation (112 nodes).Initially, material 1 is maintained with all nodes restricted and the free face of material 2 (elastic material) is stretched from coordinate 1m z = through coordinate 3.4m z = using 120 equally spaced steps.Figure 29 shows the achieved configuration with the Cauchy elastic stress distribution 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)m of material 1 and the elastic Cauchy stress at point (0.167; 0.167; 0.833)m of material 2. These coordinates correspond to the initial configuration.
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 2 1m through 2 0.31m which explains the large difference between the initial viscous stress and the final elastic stress of material 1.The adopted tolerance is As one can see by this example, the formulation is capable to simulate viscoelastic (Kelvin-Voigt-like model) structural elements including dynamics.

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.
which means that vol E is the Lagrangian strain direction corresponding to the hydrostatic stress in the Cauchy space and the hydrostatic second Piola-Kirchhoff stress intensity.The identification of the Lagrangian hydrostatic stress vol elas S and the corresponding strain direction vol E 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 ( 1 I ) of the isochoric part of the Cauchy-Green stretch tensor C as:

Fig. 1 :
Fig. 1: Deformation function described by prismatic elementsThe position based kinematics description is given in a brief way as more details can be seen in referencesSanches and Coda (2013),Coda (2015),Coda HB and Paccola (2011) andPascon and Coda (2013) for instance.In Fig.1one observes that the initial configuration mapping describes the initial coordinates of continuous points i x as a function of 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): (2005) treats the fluid as incompressible we adopted a high value for the bulk modulus ( 9 2.15x10 K =

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

Fig. 5 -
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 -
Fig. 6 -Displacement fields of dam-rupture model using 392 cubic elements at different times: (a) * 0.896 t = , (b) * 2.391 t = , wave is detected by exhibiting the first 20 time steps of this analysis, being the first static.Remembering that the pressure wave velocity

Fig. 12
Fig.12shows the vertical displacement of the upper right and left points for the two discretizations, without considering surface tension.Fig.13shows the effect of the water surface tension (

.
The application of this surface tension is made as explained in the previous example and its corresponding nodal force value is et al. (2011), for instance.The adopted time interval is 10ms .Two types of analyses are preformed, without any additional damping (for which the bubble vibrates indefinitely) and other introducing an artificial mass proportional damping ( 3 1kg.(ms)/mc = the 8x8 cubic approximation and near incompressible bulk modulus (

Fig
Fig. 15 -Bubble behavior along time 2.15GPa K = (cubic approximation) wave has a very low speed and the corresponding behavior is almost static.

Fig. 16 -
Fig. 16 -Position behavior of the analyzed bubble along time using a mesh of 4x4 cubic elements and bulk modulus 4 2.15x10 GPa K −

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

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.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).

Fig. 20 -
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

Fig. 21 -
Fig. 21 -Pressure fields at various time ( retaining the cement paste in an inverted hollow bucket, with its subsequent instantaneous release on a smooth surface.Fig.22depicts the initial configuration geometry ( 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 0.5 γ =

Figure 24 Fig. 23 -
Figure24shows the time history of the relation inf (t) / num exp final D D for the three adopted rheological models, i.e.,

=.
direction, see Fig. 25 ( ) 40kN F = and its longitudinal strain ( 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: We used 3 different values for shear viscosity:

Figure 26
Figure 26 presents the time response of the stretch in z direction and snapshots of the Cauchy viscous stresses 11 vis σ

Fig. 28 -
Fig. 28 -Indirect relaxation test: (a) Material distribution; and (b) Discretization for a quarter of the problem (Symmetry) Fig. 29 -Indirect relaxation test: (a) Initial position; and (b) Cauchy stress along time 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 σ , second Piola-Kirchhoff stress tensor S , real strain variation δε and Large strain Flory's decomposition for Lagrangian modeling of viscoleastic solids and compressive fluids Renato Takeo Kishino et al.
Latin American Journal of Solids and Structures, 2022, 19(4), In the beginning of a time step Latin American Journal of Solids and Structures, 2022, 19(4),