A mixed u–p edge-based smoothed particle finite element formulation for viscous flow simulations

A mixed u–p edge-based smoothed particle finite element formulation is proposed for computational simulations of viscous flow. In order to improve the accuracy of the standard particle finite element method, edge-based and face-based smoothing operations on the displacement gradient are proposed for 2D and 3D analyses, respectively. Consequently, spatial integration involving the smoothing operator is performed on smoothing domains. The constitutive model is based on an elasto-viscoplastic formulation allowing for simulations of viscous fluid or fluid-like solid materials. The viscous response is modeled using an overstress function. The performance of the proposed edge-based smoothed particle finite element method (ES-PFEM) is demonstrated by several numerical benchmark studies, showing an excellent agreement with analytical and reference solutions and an improved accuracy and computational efficiency in comparison with results from the standard PFEM model. Finally, a numerical application of the ES-PFEM to the computational simulation of the extrusion process during 3D-concrete-printing is discussed.


Introduction
Large deformations of solids need to be considered in numerous engineering applications, including geotechnics, fluid-structure interaction or manufacturing processes. The standard finite element method (FEM) often fails to solve problems characterized by extremely large changes of the configuration of the analyzed structures mainly due to severe mesh distortions and associated poor accuracy of the solution. Over the last decades, numerical methods have been developed to tackle these problems in a better or more elegant way. The arbitrary Lagrangian Eulerian (ALE) method [1,2] circumvents large element distortions by a combination of an Eulerian and Lagrangian description of motion. Specifically for simulations involving large motions, various particle methods such as the particle in cell (PIC) or material point method (MPM) [3], smoothed particle hydrodynamics B Günther Meschke guenther.meschke@rub.de 1 Institute for Structural Mechanics, Ruhr-University Bochum, Bochum, Germany (SPH) [4] or the discrete element method (DEM) [5] have also been developed.
In one way or another, these methods use a Lagrangian description of motion with the help of a background mesh, with kernel approximations or in a complete discrete fashion. The particle finite element method (PFEM) is also a discretization method, which uses standard Lagrangian finite elements in conjunction with a continuous re-meshing during the analysis. PFEM is a well established method that has been originally developed to solve problems in fluid dynamics [6,7] and fluid structure interaction problems [8][9][10][11]. The governing equations are solved with standard linear triangular or tetrahedral finite elements in an updated Lagrangian formulation. Hence, one advantage of PFEM is that no convective terms need to be treated in the governing equations. Beside in fluid dynamics problems, PFEM also found application in geotechnical applications, e.g. for the simulations of landslides [12][13][14][15] and excavation processes [16,17] and for solving problems in manufacturing, such as metal forming [18,19] and metal cutting processes [20]. Typically low order triangular and tetrahedral elements are used in PFEM algorithms in order to guarantee a fast and robust re-meshing of the domain. Evidently, the disadvantages of low order ele-ments, like poor accuracy and shear and volumetric locking are inherited.
More recently, the smoothed particle finite element method (SPFEM), adopting formulations underlying the class of smoothed finite element methods (SFEM) (see [21] for an overview), has been proposed. With the help of smoothening techniques, the properties of low order finite elements in PFEM models can be improved. Furthermore, depending on the specific smoothed finite element formulation, stresses are carried on the nodes or edges instead at the elemental Gauss point locations. This may reduce the error associated with the mapping of internal variables during the re-meshing procedure. The smoothed particle finite element method (SPFEM) was recently proposed by Zhang et al. [22] for geotechnical applications.
SFEM utilizes strain smoothing techniques, by smoothing the shape function derivatives over certain smoothing domains [21,23]. These smoothing domains can be defined in different ways, resulting in different SFEM formulations with different properties. The most common formulations are based on edge-based smoothing (ES-FEM) [24] or on nodebased smoothing (NS-FEM) [25]. The idea behind the SFEM is to improve the properties and accuracy of the FEM without increasing the number of degrees of freedom. SFEM effectively cures locking phenomena and overly stiff behavior of low order triangular or tetrahedral elements. Furthermore, as was shown e.g. in [23,26] SFEM improves the accuracy and convergence w.r.t. increasing degrees of freedom as well as the performance in problems with very distorted meshes. Hence, theoretically a smaller number of re-meshing steps would be necessary to guarantee a suitable discretization using an SPFEM formulation. Although the bandwidth of the stiffness matrix is larger than in standard FEM, connected with a larger computational effort for the generation of the stiffness matrix and the solution of the system of equations [26], the SFEM seems to be ideally suited to be used in a PFEM framework to cure some undesirable properties inherited from low order finite elements.
The SPFEM approach by Zhang et al. [22] for geomechanical problems is based on the node-based smoothing technique, which leads to nodal integration of terms involved in the smoothing operation and stress points located on element nodes. This model was improved in [27] by a GPU-accelerated formulation. More recently, a similar nodebased smoothing PFEM formulation was developed for fluid dynamic problems by Franci et al. [28] and was further improved for fluid structure interaction problems in [29]. The node-based smoothing formulation can effectively treat volumetric locking problems and shows upper bound properties. However, this formulation suffers from temporal instabilities in dynamic problems resulting from "zero-energy" spurious modes (see [21,30,31]). Temporal stabilization procedures may solve this issue, but may also result in a reduced accuracy of this formulation. In contrast, the edge-based smoothing formulation does not tend to be temporally unstable and stresses are carried on element edges. However, no upper bound properties are obtained with this formulation and volumetric locking is also not treated. However, the edge-based smoothing formulation improves the accuracy and convergence rates in comparison to a node-based smoothing formulation and a standard FE formulation [24]. Furthermore, the edge-based smoothing formulation is computationally less expensive compared to the node-based smoothing formulation [21]. An edge-based smoothing formulation was recently implemented into a PFEM model by Jin et al. [32].
In this paper, a smoothed finite element approach is incorporated into a PFEM formulation for solving large deformation problems. To circumvent issues related to temporal instabilities (as they would appear in node-based smoothing approaches), an edge-based approach is adopted to improve the standard PFEM formulation with low order finite elements. This edge-based smoothed PFEM formulation (ES-PFEM) results in stress points located on element edges. In order to treat problems associated with volumetric locking, a mixed u-p formulation is used. Since, due to equal order interpolation functions of the displacement and pressure field, the Ladyzhenskaya-Babuska-Brezzi (LBB) condition is violated, a direct pressure stabilization is used to address this problem. As shown in this paper, compared to the standard PFEM and the previously developed displacement SPFEM formulation, the proposed stabilized mixed u-p ES-PFEM approach considerably improves the accuracy and computational efficiency. The adopted constitutive model is based on a hyperelasticviscoplastic formulation, with the viscoplastic contribution formulated by an overstress function such as to represent a Bingham-viscoplastic response. This model is ideally suited for representing viscous fluids and fluid-like solid materials such as fresh concrete [33][34][35] or fluidized geomaterials [12,13]. As compared to regularized Bingham models [36], elasto-viscoplastic models may trace initiation and stoppage of the flow more accurately and also describe residual stresses in unloading conditions. These advantages are of substantial interest for more advanced applications, like computational simulations of additive manufacturing of concrete structures, where the material shows a rather complex response somewhere between a fluid and a solid state during the printing processes [37,38].
The structure of the remainder of the paper is as follows: In Sect. 2 the theoretical framework of the proposed ES-PFEM model is given. This section provides an overview of the constitutive model for fresh concrete, the governing equations and the weak form of the numerical model. Section 3 describes the edge-based smoothed particle finite element formulation, the finite element discretization and the algo-rithmic formulation to solve the resulting linearized system of equations. In Sect. 4 a number of numerical benchmarks are discussed to investigate the properties and features of the proposed model. Conducted numerical studies are presented for elastic benchmarks, including the (almost) incompressible Cooks membrane and an oscillating elastic beam example and viscoplastic examples, including the planar Poiseuille flow benchmark, a channel-flow test and a slump-flow test. Finally, in Sect. 4.6 a 3D-concrete-printing test is simulated numerically by means of the proposed ES-PFEM model.

Theoretical framework
The modeling strategy of the proposed formulation must be able to account for large deformations, finite strains and evolving free surfaces to adequately capture the kinematics involved. In this section an appropriate material model and the underlying strong and weak forms of the problem are introduced, which form the basis for a consistent theoretical framework at finite strains.

Constitutive model and kinematic relations
In the proposed model, the constitutive description of the material is based on a hyperelastic-viscoplastic formulation with a local multiplicative decomposition of the deformation gradient F i j = ∂ x i /∂ X j into elastic and viscoplastic parts as with the elastic and viscoplastic deformation gradients F e and F vp and with X and x denoting the coordinates in the reference and current configuration, respectively. Hence, an intermediate unstressed viscoplastic configuration is postulated, which can be obtained by a purely elastic unloading in the current configuration. Under the assumption of isotropic elastic behavior and no further internal variables, a hyperelastic potential ψ(b e ) can be defined as a function of the elastic left Cauchy-Green tensor b e , Given the multiplicative decomposition of the deformation gradient according to Eq. (1), the elastic and viscoplastic velocity gradients are obtained as Note that L vp is defined on the intermediate viscoplastic configuration. However, using Eq. (1), the spatial viscoplastic velocity gradient can be formulated as with the spatial velocity gradient l. Furthermore, the velocity gradients can be decomposed into symmetric and skew symmetric tensors as Based on the kinematic relations given in Eqs.
(1)-(6), the constitutive model is outlined in more detail in the following sections using a formulation in terms of the Kirchhoff stresses τ . The Kirchhoff stresses are related to the true (Cauchy) stresses σ via τ = σ J , with J denoting the determinant of the deformation gradient.

Viscoplastic model
An elasto-viscoplastic model is proposed to represent the behavior of Bingham-like fluid models, allowing for simulations of fluid or fluid-like solid materials with a single constitutive law. The adopted viscoplastic formulation is based on a representation of the plastic multiplier using an overstress function [39]. Hence, the algorithmic scheme of the visoplastic model is similar to standard rate independent plasticity models. The material need to satisfy the von Mises yield function with the yield stress τ 0 and the equivalent stress τ eq = 1 2 dev[τ ] : dev[τ ]. The rate of the left Cauchy-Green tensor (Eq. (2)) is given as [40] with the Lie derivative of the left Cauchy-Green tensor L v b e . Assuming material isotropy, and using Eq. (8), the plastic dissipation can be written in the following form With the assumption of the plastic spin W vp = 0, the Lie derivative is given as with d vp specified by the flow rule [41] d vp =γ m, (11) withγ denoting the plastic multiplier and m the direction of the plastic flow. Adopting an associative flow rule following from the principle of maximum plastic dissipation, m is given as The plastic multiplierγ needs to satisfy the Kuhn-Tucker loading/unloading conditions, defined for classical rateindependent plasticity as [42] along with the consistency conditioṅ The viscoplastic flow of the material is defined by expressing the plastic multiplier using an overstress function as [39] γ with the viscosity parameter μ and the McCauley brackets "< . >" defined as The viscoplastic Bingham flow model is obtained by setting the overstress function equal to the yield function Eq. (7), which yields

Hyperelastic potential
Considering isotropy of the material, the viscoplastic model is formulated in principal stretches according to [40]. This formulation has the advantage that all features of classical infinitesimal plasticity models are retained and standard return mapping algorithms from small strain plasticity can be used without any modification. To this end, we make use of the spectral decomposition of the elastic left Cauchy-Green tensor with the eigenstretches λ e A and the eigenvectors n A (with A = 1, 2, 3). The principal logarithmic elastic eigenstretches are then defined as ε e A = ln(λ e A ). The hyperelastic potential using principal logarithmic eigenstretches is given as with G as the shear modulus, K the bulk modulus and J e denoting the determinant of the elastic deformation gradient with ln(J e ) = ε e 1 + ε e 2 + ε e 3 . The principal Kirchhoff stresses can be obtained from Since the Kirchhoff stresses share the same eigenvectors as the elastic left Cauchy-Green tensor, The spatial constitutive tangent moduli are given as with n ABC D = n A ⊗ n B ⊗ n C ⊗ n D . The constitutive tensor is obtained from a lengthy derivation using the rate forṁ S = C :Ė with the Second Piola Kirchhoff stresses S, the material constitutive tensor C and the Green-Lagrange strain tensor E. The full derivation and a treatment for the case of equal eigenvalues based on L'Hospital rule can be found in [43][44][45].

Governing equations
The local form of the balance of momentum on the domain Ω in a time interval [0, T ] is given as where ρ denotes the density of the material, ρb i the external body forces acting on the domain and a i denotes the acceleration. The displacements u i are prescribed on the Dirichlet boundary Γ u as with the prescribed displacements u * i . The tractions t * i are prescribed on the Neumann boundary Γ t as with the normal vector n j . In case of plastic flow of the material, the flow rule (Eqs. (11) and (12)) leads to a purely deviatoric flow direction. Therefore, J = J e and plastic incompressibility arises, when the yield stress of the material is exceeded. It is well-known that a standard displacement finite element formulation produces poor results in case of incompressibility [46]. To avoid problems associated with incompressible material behavior, a mixed displacement pressure formulation is adopted. Hence, the Cauchy stresses are split into a deviatoric and a volumetric part, where I i j = δ i j , δ i j being the Kronecker delta. The stressesσ are obtained from the uncoupled hyperelastic potential, Eq. (19). Hence, the deviatoric portion of σ are directly evaluated from the elastic left Cauchy-Green tensor b e of an integration point. However, the pressure p is obtained from the additional relation [14] which can be further re-arranged using the constitutive model Eq. (19) in the form Equations (22) and (27) constitute the balance equations for the proposed smoothed finite element formulation.

Weak form
The weak form of the balance of momentum is obtained from the typical procedure by multiplying Eq. 22 with an arbitrary test function δu and integration over the current configuration as [44,47] where δ D i j = 1 2 ∂δu i /∂ x j + ∂δu j /∂ x i are the virtual deformation rates and Ω t is the volume in the current configuration. Equation (28) is already reformulated by applying the divergence theorem to the part involving σ i j and by using the traction boundary conditions Eq. (24).
Using the arbitrary test function δp, the weak form of Eq. (27) is obtained accordingly as where Ω 0 is the volume in the reference configuration. The choice of taking an integral over the reference or current configuration is arbitrary and should result in the same approximated pressures. However, integration over the reference configuration was chosen, since linearization of the integral as given in Eq. (29) is much simpler. Integrals can be transformed to another configuration by the standard relation

Numerical formulation
In this section, an overview of the numerical formulation is presented. First, a short introduction in the particle finite element method is given. Then the edge-based gradient smoothing technique is outlined and the pressure stabilized smoothed finite element approximation is described. Finally, the incremental formulation and the solution algorithm are briefly summarized.

The smoothed particle finite element method
In the original particle finite element method (PFEM), the domain is typically discretized by standard low order triangular or tetrahedral finite elements [6,7] in an updated Lagrangian description. The finite element discretization is frequently re-meshed, when the mesh becomes too distorted. This re-meshing procedure is based on a Delaunay triangulation in conjunction with the alpha shape method for detecting the free surface of the domain [48]. However, in the context of this work, two features of this original method are modified to improve the method. Firstly, constrained Delaunay triangulations are used and secondly, an edge-based smoothed finite element formulation is adopted. The proposed numerical formulation is intended for applications with free surfaces remaining continuous without self-contact or excessive motions such as wave splashing. Considering these assumptions, it is found that altering the re-meshing procedure by using constrained Delaunay triangulations offers advantages in comparison to the original method in terms of free surface recognition and mass conservation, see for instance [20]. In constrained Delaunay triangulations, the free surface of the domain must be known in advance, so that it can be used as an input for the algorithm. Hence, the connectivity of the free surface of the domain must remain unaltered over the simulation time, or it must be adapted manually. Refining, updating or keeping track of the free surface is straightforward in 2D. In 3D however, due to the higher degree of geometrical freedom, guaranteeing for a well-shaped surface triangulation is very difficult. Therefore, when necessary, the free surface is recovered from an additional Delaunay triangulation, improved by a deletion of all elements that were not located inside the old discretization. Note that the initial free surface triangulation at the start of the simulation is obtained from the alpha shape method. Furthermore, by adopting a constrained Delaunay triangulation in conjunction with a mesh quality improvement algorithm, sliver and badly shaped elements are avoided, see [49]. Typically in PFEM algorithms, contact is modeled by using contact elements that are naturally formed during the re-meshing procedure. However, for a more accurate representation of the boundaries in this work, contact with the ground is modeled by means of a standard penalty method, see [6,9,20] for more details on contact treatment in PFEM simulations. However, self-contact or contact between bodies is not included in the current implementation.
A drawback of Lagrangian particle methods is that throughout the simulation procedure clogging or separation of particles might happen, which is circumvented by deletion or insertion of particles. In case of insertion of new particles, solution variables are linearly mapped onto the new particles. Even though, PFEM was developed as a particle method in the field of fluid dynamics and soon found application in problems involving solid mechanics [12,13,[16][17][18][19]. Therefore, in addition to the solution variables, also internal variables located on the Gauss points need to be mapped to the new discretization, due to changes of nodal connectivities after the re-meshing procedure. The mapping of internal variables is accomplished by a mesh-to-mesh based nearest neighbor interpolation, which is commonly used in PFEM formulations [13,14,19,20]. With this algorithm, diffusion of internal variables during re-meshing is minimized.
The algorithmic steps involved in the smoothed particle finite element method (SPFEM) are summarized for one time step (n → n + 1) as follows: 1. Define the point cloud with the free surface and initialize the calculation of the current time step t n . 2. Generate a new finite element discretization using the constrained Delaunay triangulations of the point cloud and project all variables from the old to the new discretization. 3. Solve the system of equations in the updated configuration at t n+1 = t n + Δt. 4. Update the particle positions based on the calculated displacement increment from the previous step. Go back to step 1 and repeat this procedure for the next time step, if t n < t max .

Edge-based gradient smoothing
In contrast to elemental integration, as employed in standard FEM, in edge-based smoothing techniques the integration of terms involving the smoothing operator are performed on cells associated with each element edge. These smoothing cells are constructed from the finite element mesh by using the mid points of the elements and nodes as shown in Fig. 1.
In this work an edge-based smoothing approach is adopted for the 2D case. The adopted 3D extension of this approach is face-based smoothing. For convenience, in the following the term "edge-based smoothing" is used in all derivations, knowingly that the 3D extension is a face-based approach. However, edge-based smoothing can also be used in 3D [26], leading to even better properties than the face-based alternative. The edge-based smoothing operation is applied to the gradient of the displacement field on each smoothing domain Ω k of edge k as with the displacement gradient u i, j = ∂u i /∂ x j and the smoothing function W k . The smoothing function W k must satisfy the partition of unity. A widely used smoothing function is the Heaviside smoothing function defined as with V k denoting the area in 2D or volume in 3D of the smoothing domain Ω k . The smoothed displacement gradient is reformulated by applying the divergence theorem to Eq. (30) and using Eq. (31) as with the normal vector n k on the boundary Γ k of the smoothing domain Ω k .

Finite element discretization
The governing equations are discretized by low order triangular (in 2D) and tetrahedral finite elements (in 3D) in an updated Lagrangian formulation. Hence, the solution variables and test functions are discretized by finite element shape functions as [46,50], with the shape functions for the displacement and pressure fields N u and N p , the nodal values of the displacements u and of the corresponding test function δū, and the nodal values of the pressurep and of the test function of the pressure field δp.
Using the finite element discretization introduced in Eq. (33), the smoothed displacement gradient Eq. (32) can be formulated as where N k refers to all nodes connected to edge k.b k I j are the smoothed derivatives of the shape functions, The smoothed derivatives of the shape functions Eq. (35) can be evaluated using Gauss integration over all boundary segments of Γ k . Alternatively, when triangular or tetrahedral elements with linear shape functions are used, the values of b k I j are constant and can be obtained simply by a volume weighted average of the associated elemental derivatives of the shape functions, see [21].
The finite element discretization Eq. (33) together with edge-based smoothing of the displacement gradient Eq. (35) are applied to the weak form Eqs. (28) and (29) to yield the discretized system of equations as with the nodal values of accelerationā. The matrices and vectors in Eqs. (36) and (37) are obtained from assembling the element-and edge-based contributions as with the elemental volume Ω e and the smoothed gradient operatorB k . The mass matrix M e uu is used in lumped form. The term F e p is evaluated at the mid point of edge k. Note that all integrals are evaluated in the updated configuration. The smoothed gradient operatorB k depends on the number of nodes N k connected to edge k and is given as For the 2D case,B k I is given for a node I connected to edge k as The size ofB k varies, for internal and free surface edges, and for 2D and 3D problems, respectively, since the number of nodes N k connected to edge k changes.

Stabilization
Given the residuals Eqs. (36) and (37), using equal order linear finite element shape functions for both the displacement and the pressure field, will violate the Ladyzhenskaya-Babuska-Brezzi (LBB) condition. In this work, a direct pressure stabilization technique based on the polynomial pressure projection (PPP) is chosen to stabilize the mixed formulation [51]. The stabilization technique is based on a projection of the C 0 approximated pressure field onto a pressure field approximated with a one order lower approximation, which is consistent with the discontinuous stress field. Hence, the inconsistencies between an equal order interpolation for the displacement and pressure field are resolved. This projection operator is added to the balance of mass, leading to where α stab denotes a parameter for tuning the level of stabilization andp denotes the discontinuous approximation of p onto a polynomial space of one order lower than the approximation for p. Using this stabilization technique, no mesh dependent stabilization parameter needs to be introduced and no higher-order derivatives need to be calculated. The stabilization can be computed for each element from the definition of the projection operator as Therefore, Eq. (37) is modified as with the stability matrix K pp,stab in the current configuration and re-written by using Eq. (46) as whereN p are the discontinuous polynomials given in 2D as

Algorithmic formulation and stress update procedure
To solve the constitutive model as outlined in Sect. 2.1, an incremental kinematic procedure must be introduced. Let F n be the deformation gradient at time t n and let F n+1 be the deformation gradient at the updated time t n+1 , the incremental deformation gradient is defined by a multiplicative decomposition of the deformation gradient as with Δu n+1 as the incremental displacement and ∇ n = ∂/∂ x n . The temporal discretization of the viscoplastic part is based on an exponential integrator using the viscoplastic velocity gradient Eq. (4) as After straightforward manipulations, the updated elastic left Cauchy-Green tensor is obtained as with the trial elastic left Cauchy-Green tensor and its representation in principal stretches given as Note that b e n+1 and b e,trial n+1 share the same eigenvectors n trial A by definition, see for instance [40]. Furthermore, by applying the tensor logarithm to Eqs. (51) and (52) the (trial) principal logarithmic stretches are given as ε e A = ln(λ e A ), which are updated according to where m A,n+1 is the direction of the plastic flow (see Eq. 12) in principal axes A = 1 . . . 3. Hence, the combination of an exponential integrator and a hyperelastic model based on principal logarithmic stretches leads to an additive decomposition into elastic and viscoplastic parts of the principal logarithmic strains and, therefore, a standard small strain return map algorithm can be used when the yield function is violated, for details we refer to [40][41][42].
The updated elastic left Cauchy-Green tensor and Kirchhoff stresses are then defined as The updated Cauchy stresses are further obtained with σ n+1 = τ n+1 /J . Consequently, the only internal variable defining the current stress state is the elastic left Cauchy-Green tensor b e , which, as a reminder, is defined on the element edges in 2D or faces in 3D due to the adopted smoothed FE formulation. As outlined in Sect. 3.1, the internal variables are projected by a mesh-to-mesh based nearest neighbor interpolation during re-meshing. However, it was found that during remeshing and after projection of internal variables the new residual R p in Eq. (47) may be perturbed away from equilibrium, which would slow down the convergence of the solution procedure, and more iterations would be necessary. Accordingly, the projection operation is relaxed with respect to the determinant of the deformation gradient J by replacing the projection of J by an incremental superconvergent patch recovery projection to yield a smoother residual R p . Therefore, the determinant of the smoothed incremental deformation gradient j k n+1 = det( f k n+1 ) of an edge k is projected on the node l, using the volume average as follows: with j l n+1 being the relative determinant of the deformation gradient of node l, N l are all edges connected to node l and h is 2 or 3 in 2D or 3D, respectively. During re-meshing, the isochoric elastic left Cauchy-Green tensor is still obtained from the original projection algorithm, while the volumetric part of the elastic left Cauchy-Green tensor is recovered from the nodal determinant of the deformation gradient as given in Eq. (56). This procedure considerably improves the convergence properties and renders the solution process more robust, while the introduced diffusion effect remains small. A numerical example demonstrating the effect of this modification is discussed in Sect. 4.2.

Solution algorithm
Equations (36) and (47) define the final system of equations, which must be solved iteratively. Note that discretization in time of the kinematic variables is performed by the implicit α-Bossak scheme [52]. In the following, the solution algorithm is summarized for one time step (n → n +1) to obtain the results in the updated configuration at t n+1 . The iteration counter is denoted as i defining the ith iteration. The solution procedure is described here, for simplicity, for the 2D version. The implementation has been performed in 2D and 3D, with the extension to 3D being straightforward. It can be partitioned into 7 algorithmic steps: (1) Initialization of variables: Setx i n+1 =x n ; Δū i n+1 = 0; p i n+1 =p n ; b e,i n+1 = b e n (2) Calculate the residuals R i u,n+1 and R i p,n+1 according to (36) and (47).
(3) Calculate the tangent matrices as where the sub-matrices are given as [44] with α and β denoting time integration parameters from the α-Bossak scheme [52], m = [1, 1, 0] T in 2D and the elemental geometrical stiffness matrix K e uu,geo and the material stiffness matrix K e uu,mat are defined as with the algorithmic constitutive modulus of the mixed formulationc k = c k dev,evp + p(I ⊗ I − 2I), where I i jkl = (δ ik δ jl + δ il δ k j )/2 is the fourth order symmetric identity tensor. The constitutive tensor c k dev,evp is simply the deviatoric version of Eq. (21) in conjunction with the consistent small strain elastic-viscoplastic consitutive tensor [44,45]. The matrixG k is (in 2D) given for a node I connected to edge k as

(5) Update the solution variables
and update the velocities and accelerations according to the α-Bossak scheme [52]. (6) Update the elastic left Cauchy-Green tensor and the stresses according to Sect. 3.5 and perform a return map algorithm, when the yield condition Eq. (13) is violated.
(7) Calculate the residuals using the updated variables by repeating step (2) and check for convergence: When R i+1 u,n+1 / R u,n < tol and R i+1 p,n+1 / R p,n < tol, the solution has converged and the time step calculation is finished. Set n = n + 1 and repeat the calculation for the next time step by starting at step (1). If the solution has not converged, update the iteration counter by setting i = i + 1 and repeat steps (3)- (7). In this work the value of the error tolerance tol was set to 10 −10 .
Due to the adopted edge-based strain smoothing technique, the bandwidth of the stiffness matrix is a few times larger than in standard FEM using the same discretization and therefore the overall solution algorithm is computationally more expensive, see [21]. Hence, it was found that calculating the tangent matrices only once at the beginning of each time step using the values from the previous time step was the most efficient way to solve the system of equations in the studied examples in this work. Due to generally small time steps employed in PFEM simulations, the solution typically converged within a few iterations. Further improvements and modifications may be tested in the future, including staggered schemes [10,11,53] or explicit algorithms [54]. Also note that contributions from the contact definition are omitted here for simplicity. As already pointed out in Sect. 3.1, a standard penalty method is adopted for contact with the ground surface.

Numerical benchmark studies
In this section different numerical studies are conducted to evaluate the performance of the proposed edge-based smoothed PFEM (ES-PFEM) formulation. In the first benchmark study the Cook's membrane is used to demonstrate the properties of the proposed numerical model in comparison to other numerical formulations (e.g. a pure displacement FE formulation) for a quasi-incompressible problem. The second study shows the performance of the proposed numerical model in a dynamic problem, using an elastic dynamic cantilever beam example. The first two studies reveal the superior properties of the proposed edge-based smoothed FE model compared to a standard FE formulation. In the third example the planar Poiseuille flow is used to verify the implemented constitutive law to simulate the flow of a Bingham fluid. The next two studies are the channel flow test and the slump flow test, which are concrete flow related benchmarks. The channel flow test is studied in 2D, while the numerical simulation of the slump flow test is carried out in 3D. The final example is concerned with the modeling of the extrusion process of a 3D-concrete-printing scenario. are investigated and compared in the presented studies. The standard FEM models are based on the same theoretical framework as the numerical features outlined in the previous sections.

Cook's membrane
In this numerical study the properties of the proposed smoothed finite element formulation are analyzed and compared with different finite element formulations. This analysis shows the superior performance of the proposed model for solving (quasi-)incompressible problems. The Cook's membrane is a well-known benchmark for testing finite element formulations in 2D [55]. In accordance with the original benchmark, a 2D plain strain model was set up in a geometrically linear setting. Therefore, all nonlinear contributions in the balance equations were omitted and a standard small strain linear elastic material model was used, considering a quasi-incompressible material with ν = 0.49999 and E = 250. The geometry and remaining parameters are presented in Fig. 2. The structure is discretized by three edge-based smoothed finite element formulations (ES-FEM-USP, ES-FEM-UP, ES-FEM-U) and two standard FE formulations (FEM-UP, FEM-U) for different levels of discretization.
In Fig. 3 the displacement at the top right corner of the investigated membrane is plotted as a function of the number of degrees of freedom for the different finite element formulations and discretizations. For comparison, the reference solution taken from [56] was included as a dashed red line. As can be observed, the pure edge-based smoothed displacement formulation (ES-FEM-U) tends to show volumetric locking and the results are only slightly better than the standard FE displacement formulation (FEM-U). This locking phenomenon can be effectively treated by the incorporation of a coupled displacement pressure formulation (ES-FEM-UP), which increases the accuracy of the results dramatically. Furthermore, by additionally stabilizing the pressure (ES-FEM-USP, α stab = 0.5), the accuracy is even further improved. The final results are also in a good agreement with the reference solution [56] and only vary by less than 1%. Furthermore, the results obtained with edge-based smoothing are much more accurate in comparison to the results obtained with the standard FE formulations.

Elastic dynamic cantilever beam
This study serves as a benchmark example for solving dynamic problems. In this analysis the proposed model (ES-FEM-UP) is also compared with standard FE formulations without edge-based smoothing.
The study in focus is an undamped oscillating short elastic cantilever beam subjected to a shear stress on the right edge in 2D plain strain, see Fig. 4 for more details about the geometry. The material parameters are E = 4 GPa, ν = 0.32 and ρ = 1450 kg/m 3 . This problem was studied in [57] and compared to an analytical solution of a short beam. Furthermore, this study was already carried out for a PFEM formulation in [8]. The simulation was studied with a time step size Δt = 0.001 s and a maximum time of t max = 1 s. Different finite element discretizations (5 × 17, 9 × 33, 17 × 65 and 33 × 129 number of nodes in height and length direction) were analyzed to assess convergence of the results. Furthermore, the displacement of the mid-point on the right edge over time is given in Fig. 6 for different discretizations and the proposed velocity-pressure edge-based smoothed finite element formulation (ES-FEM-UP). The results predict a good agreement with the analytical solution provided in [57], which is also given in Fig. 6. The same results for an equivalent standard FEM formulation (FEM-UP) are depicted in Fig. 7. In comparison to Fig. 6 it can be observed that the proposed edge-based smoothed formulation gives more accurate results, especially for very coarse discretizations.
In order to analyze the influence of the modified projection algorithm of J according to Eq. (56), the incremental This can also be assessed from Fig. 9, where the maximum displacement is given for different finite element discretizations and formulations (including pure displacement formulations). Results obtained from a coupled velocitypressure formulation in comparison to a pure displacement formulation also give higher accuracy. However, this effect is stronger for the standard FE formulations (FEM-U and FEM-UP). This indicates that shear locking problems are Finally, based on computations on a standard personal computer, the relative error of the maximum displacement at the mid-point on the right edge is plotted over the accumulated simulation time for different finite element formulations in Fig. 10. The displacement error was calculated as e u = (u max − u max, f d )/u max, f d with u max, f d as the maximum displacement obtained from the finest discretization of the respective finite element formulation. In accordance with observations reported in [24], the CPU-time versus error diagrams demonstrate the superior computational efficiency of the edge-based smoothing approaches as compared with the standard FE formulations.

Planar Poiseuille flow
The Poiseuille flow is a well-known benchmark problem in fluid dynamics conducted as a verification test for the numerical model including the constitutive model presented in Sect. 2.1. The geometry of the planar Poiseuille flow problem is depicted in Fig. 11. No-slip boundary conditions are applied along the walls and the Poiseuille flow of a Bingham fluid for three values of the yield stress was studied. According to [58] the analytical solution of the planar Poiseuille flow of a Bingham fluid for the velocity in x-direction over with the dynamic yield stress τ 0 , the plastic viscosity μ, the channel height h, the body force F x and the channel height z t = τ 0 /F x at which the material yields. The model and material parameters were specified as μ = 10 Pa s, τ 0 = 20, 50, 80 Pa, ρ = 1000 kg/m 3 , F x = 20,000 N/m 3 and h = 0.2 m. The elastic parameters were chosen large enough, so that elastic strains became negligible. The simulation was carried out by neglecting inertia effects to reach the steady state solution in one time step. The computed velocities in x-direction over the channel height z are plotted together with the analytical solutions for the three levels of yield stress in Fig. 12. The characteristic plug flow at the channel center, where the yield stress τ 0 of the fluid is not exceeded, as well as the zero shear flow of the material at the channel walls are well replicated and an overall perfect match between analytical and numerical results is observed.

Channel flow test
In the following numerical study a benchmark test related to the flow of fresh concrete is simulated using the proposed ES-PFEM model. The conducted numerical analysis is based on the channel flow test, which is a suitable laboratory test to asses the yield stress of self compacting concrete (SCC). The experiment is characterized by the following procedure [34]: Over 30 s 6 l of fresh concrete are carefully poured into one end of a rectangular channel (l = 1.2 m, w = 0.2 m, h = 0.15 m). After this pouring phase, a creep flow of the material develops. When the material finally comes to rest, the resulting shape of the fresh concrete allows the yield stress τ 0 of the material to be estimated. This laboratory test is also used as a numerical benchmark problem for calibrating numerical models. In [35] this test was performed in a comparative evaluation of different numerical models for a fresh concrete with the following material properties: μ = 50 Pa s, τ 0 = 50 Pa and ρ = 2300 kg/m 3 . The same numerical test was re-analyzed in this example with E = 0.1 MPa and ν = 0.3. The results obtained for the final shape of the concrete were compared with the analytical solution provided in [34]. The numerical analysis was performed in 2D assuming plain strain conditions with different mesh discretizations (1300, 4800 and 7300 elements at the end of the simulation), a time step size of Δt = 2 × 10 −4 s and a maximum time step of t = 200 s. The stabilization parameter was chosen as α stab = 200 to yield a stable pressure field. Note that the effective stabilization is approximately the same as in the Cook's membrane problem due to scaling of the stabilization term in (45) with the shear modulus G. The same stabilization parameter was also used in the remaining studies. The pouring process of the material was simulated by an inlet flow such that new elements were continuously pushed into the channel, see Fig. 13.
The simulated shapes at different time steps (t = 5 s, t = 15 s, t = 30 s and t = 200 s) obtained for the finest finite element discretization are given in Fig. 13. In Fig. 14 the velocity magnitude is depicted for the time steps t = 5 s, t = 15 s, t = 30 s and t = 200 s. As can be observed, the velocity decreases to a very low number in the order of 10 −5 m/s at t = 200 s, indicating an extremely slow creep flow at t = 200 s. However, since the final shape of the material is only influenced slightly by such a small velocity, we assume stoppage of the flow at this point. Similar observations can be made from Fig. 15, which shows the flow distance of the concrete material over time for different finite element discretizations. As can be seen, the flow almost stops at around 80 s. Furthermore, it is noted that all three discretizations predict almost identical results. The final flow distance at 200 s is 0.70532 m for the finest discretization. This value varies by 1.02% from the average value found in [35] for different numerical methods. Finally, a perfect agreement of the final shape for the finest discretization in comparison to the analytical solution, provided by [34], is found, see Fig. 16.

Slump flow test
Similar to the channel flow test in the previous Sect. 4.4, the slump flow test is also a practice-oriented test to assess the workability of fresh concrete. A conical frustum, see Fig. 17 is placed on a flat surface, filled with the fresh concrete and then lifted up, so that the material can flow in all directions.   [34] The final spread of the material gives a measure for the workability of the tested fresh concrete [33]. The slump flow test was already studied in a PFEM framework in [59,60] and SPH framework [61]. However, in this analysis the same concrete specifications as provided by the example in [35] and used in the previous Sect. 4.4 (μ = 50 Pa s, τ 0 = 50 Pa, ρ = 2300 kg/m 3 , E = 0.1 MPa and ν = 0.3) was analyzed. The final shape of the slump obtained from the proposed ES-PFEM model was compared with the analytical solution given in [33].
This numerical analysis was conducted in a 3D framework by simulating one quarter of the cone, exploiting the symmetry of the problem. A no slip condition was applied to the ground surface. The analysis was performed with three different finite element discretizations (10,000, 24,000 and 79,000 elements at the end of the simulation). The time step size was Δt = 5 × 10 −4 s with a maximum time step of t = 60 s.
The shape of the slump at four stages of the setting (after t = 0.1 s, t = 0.2 s, t = 1 s and t = 60 s) is shown in Fig. 18 for the finest mesh discretization. Figure 19 shows contour plots of the magnitude of the velocity at three stages after setting (t = 0.2 s, t = 0.5 s, and t = 60 s). While in the initial stage, the maximum velocity (in red color) is distributed at the top of the conetype shape, the flow distribution changes towards the edges of the slump as it continues to spread radially. Similar as for the channel flow test, a very slow creep flow with a velocity in the order of 10 −4 m/s is observed at 60 s. Hence, the final spread of the material is only influenced slightly after t = 60 s, as is also confirmed by Fig. 20, which contains the radial spread of the fresh concrete over time for three different mesh discretizations. The final shape of the cross section of the slump at t = 60 s is compared with the analytical solution [33] in Fig. 21. The computed diameter of the proposed model varies by 6.01% in comparison to the conducted numerical analyses in [35] and the average value found therein.

Numerical simulation of 3D-concrete printing
In this section, the ES-PFEM model is applied to 3D computational simulations of the extrusion process of a concrete printing experiment. In concrete printing, the fresh, flowable concrete is extruded layer-wise to form structural components without any use of additional scaffolding support during casting [62,63]. It was shown in [63][64][65] that the material deformation of fresh concrete during extrusion can be modeled with sufficient accuracy using viscoplastic consti-  [33] tutive models. The conducted virtual experiment is based on a one layer test documented in [64]. The material was printed with a round extrusion nozzle with a diameter of 5.6 cm, a printing speed of 0.05 m/s, a distance between the extrusion nozzle and the ground of 2.4 cm and a printed layer length of 50 cm, see Fig. 22. The material flow inside the extrusion nozzle was not included in the simulation. Hence, the extrusion flow of the fresh concrete was applied directly at the outlet of the extrusion nozzle with a constant inlet velocity of 0.03447 m/s. To retain the full Lagrangian description of the model, the inlet flow was realized by continuously pushing new elements into the domain. Furthermore, a no-slip condition was applied to the interface between the fresh concrete and the ground surface. Based on a previous numerical-experimental analysis [64], typical material properties of 3D-printed fresh concrete were adopted as follows: μ = 5 Pas, τ 0 = 306 Pa, ρ = 2078 kg/m 3 , E = 0.1 MPa and ν = 0.3. The time step size was specified as Δt = 2 · 10 −4 s. The material at the extrusion nozzle was approximated with three different finite element discretizations using a node spacing of approximately 11, 15 and 19 number of nodes along the diameter, which resulted in approximately 16,000, 80,000 and 184,000 finite elements after t = 10 s.
The printed shapes obtained from the finest mesh discretization are shown in Fig. 22 for three stages of the extrusion process (t = 2 s, t = 6 s and t = 10 s). As can be observed, already after few seconds of printing, a steady state printed cross section is obtained. The printed cross sections resulting from analyses using three mesh discretizations are shown in Fig. 23. The shapes obtained for the different discretizations almost coincide. This implies that already a coarse spatial discretization is able to predict the expected cross section with sufficient accuracy. The model results obtained for the finest mesh discretization are validated by comparing the printed cross-section with measurements from laboratory tests and with results from numerical simulations using the regularized Bingham model, both taken from [64]. Figure 24 shows the shapes of the printed concrete after 10 s. Almost no difference is found in comparison with the regularized Bingham model, and only marginal differences are observed when comparing the numerical predictions and the experimental data.
In Fig. 25 the spatial distribution of the magnitude of the velocity (top) and the pressure (bottom) after t = 10 s is illustrated in Fig. 22 along a longitudinal section and the cross section, respectively. Evidently, the maximum velocity appears directly below the extrusion nozzle, while quickly reducing almost to zero when the material touches the ground surface. As already mentioned, after a few seconds of printing  Fig. 24 3D simulation of concrete printing: simulated cross section after t = 10 s. Comparison between the proposed ES-PFEM model and the regularized Bingham model and experiments [64] time, a steady state is obtained. The figure also shows that the extrusion flow causes a pressure distribution, with the maximum pressure spreading from the rear part of the extrusion nozzle to the ground surface. This pressure distribution along the bottom face of the printed layer may significantly influence the stability and the deformations of the printed structure during multilayer printing. Figure 26 shows the distribution of the equivalent stress after 10 s obtained from the proposed ES-PFEM formulation (top) and the standard PFEM formulation (bottom) along the longitudinal section and the cross section, respectively. The equivalent stress exceeds the yield stress only directly below the extrusion nozzle. Hence, immediately after extrusion of the material, the material response is governed by the regime below the yield stress. In contrast to regularized Bingham models often used for the flow of fresh concrete, the proposed elastic-viscoplastic model is able to capture the stresses in this regime. This can be advantageous during multilayer printing scenarios, when elastic deformations must be captured, which may result in residual stress states during printing. By comparing the equivalent stresses with the standard PFEM model in Fig. 26, a more continuous distribution of the stresses is observed for the proposed ES-PFEM formulation. However, stresses in the ES-PFEM formulation are projected from the faces to the elements for a more convenient post-processing, which may also contribute to some extent to the slightly smoother stress state. Nevertheless, measuring the influence of re-meshing and stress projections in PFEM is not a trivial task. Based on the visual inspection of Fig. 26, when adopting the smoothed FE formulation, a more smooth stress state is observed which seems to be less influenced by the re-meshing and stress projection procedure.  The numerical formulation is inspired by Smoothed Finite Element Methods (SFEM) in order to improve the performance of low order finite elements typically employed in PFEM algorithms. An edge-based smoothing operation was applied to the displacement gradient. This technique results in stresses being carried on edges, instead of Gauss points as in standard PFEM formulations. Hence, a more continuous stress field over the element edges was enforced resulting in higher accuracy of the numerical analysis results in comparison to standard low order elements. The adopted constitutive model was based on a hyperelasticviscoplastic consitutive law as it allows for the modeling of viscous fluids and fluid-like solids. The hyperelastic potential was defined using principal logarithmic stretches and the viscoplastic contribution was modeled by a Perzyna-type overstress function. Due to the purely deviatoric flow rule and the associated plastic incompressibility, a mixed u-p formulation was used to archive high accuracy with the given model.

Conclusions
The discussed examples in this work demonstrated the capabilities of the numerical model and the superior properties of the proposed ES-PFEM formulation. In the first study the Cook's membrane was analyzed for the (quasi-)incompressible case, by comparing convergence and accuracy of different numerical formulations. This example revealed that the proposed model can efficiently avoid volumetric locking and is therefore suitable for the proposed viscoplastic formulation with a purely deviatoric flow rule. In the second analysis an oscillating elastic cantilever beam was taken as an example to compare the original PFEM model with the proposed ES-PFEM model. This study showed a higher accuracy and computational efficiency of the proposed ES-PFEM model in comparison to standard PFEM models. Hence, a higher accuracy was obtained with less CPU-time consumption. Moreover, it was shown that the incremental superconvergent patch recovery of J according to Eq. (56) (used as a projection procedure of J during re-meshing in PFEM) resulted in negligible diffusion and increased accuracy of the formulation. The first two studies indicated that the inherited properties of low order finite elements in PFEM can be efficiently improved by the proposed edge-based gradient smoothing technique. The third study served as a verification test for the elasto-viscoplastic model to be able to replicate the behavior of a Bingham fluid model. The analytical solution of the planar Poiseuille flow test was compared with the proposed model, and a perfect match was found. Finally, the fourth and fifth numerical benchmark studies were concerned with simulations of the channel flow test in a 2D framework and of the slump flow test in a 3D setting. Both benchmark simulations showed good correlations with the analytical solutions, which demonstrates the ability of the proposed numerical model to adequately simulate the flow behavior of fresh concrete. Interestingly, almost no influence of the mesh density on the predicted flow behavior was found. In the final computational study the extrusion process of a concrete printing experiment was analyzed by means of the proposed ES-PFEM model in a 3D setting. The calculated shapes of the extruded concrete layer were compared to laboratory tests from the literature. A good agreement was found even for the coarsest discretization. Due to the elastic-viscoplastic description of fresh concrete, stresses in the regime below the yield stress can be captured. Hence, this model has the ability to track elastic deformations and resid-ual stresses during multilayer printing processes. For realistic simulations of additive manufacturing of concrete at a larger time scales, however, a more elaborated model is required to account for structural build-up, plastic shrinkage and creep which are relevant physical processes in printed concrete in the deposition scale of 3D-concrete-printing [63].
In conclusion, it was shown that the proposed ES-PFEM formulation is superior over the standard PFEM with respect to accuracy and computational efficiency. The proposed formulation also showed minimal diffusion due to projection operations of J during re-meshing and is potentially less sensitive to mesh distortions [21], theoretically leading to a smaller number of re-meshing steps. In addition, the proposed implementation requires only minor changes in an existing PFEM code. Despite the superior properties of the proposed ES-PFEM model in comparison to standard PFEM formulations, smoothed finite element formulations usually require a larger amount of computational effort due to a larger bandwidth of the stiffness matrix using the same discretization. For this purpose it might be desirable to rely on a fully explicit algorithm in the future, to avoid the calculation of the stiffness matrix. It is noted that the projection of internal variables during re-meshing is still necessary for the proposed formulation. In this context, further investigations are necessary and it might also be advantageous to further investigate different smoothing techniques, like imbricate [66,67] or node-based smoothing approaches, to avoid re-meshing related projections of internal variables [22,28], which are associated with individual unique properties. However, node-based SFEM models require a temporal stabilization potentially leading to reduced accuracy. Even though stresses do not need to be projected during remeshing, there still might remain a re-meshing related error due to new nodal connectivities and, consequently, changes of the smoothing volumes. Therefore, it is still an open question, if such models would lead to an improvement of the proposed model.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.