- Email: [email protected]
A stable P1∕P1 finite element for finite strain von Mises elasto-plasticity
Accepted Manuscript A stable P1/P1 finite element for finite strain von Mises elasto-plasticity E. Feulvarch, J.-C. Roux, J.-M. Bergheau, P. Gilles
PII: DOI: Reference:
S0045-7825(17)30169-X http://dx.doi.org/10.1016/j.cma.2017.06.026 CMA 11491
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date : 24 January 2017 Revised date : 6 June 2017 Accepted date : 22 June 2017 Please cite this article as: E. Feulvarch, J.-. Roux, J.-. Bergheau, P. Gilles, A stable P1/P1 finite element for finite strain von Mises elasto-plasticity, Comput. Methods Appl. Mech. Engrg. (2017), http://dx.doi.org/10.1016/j.cma.2017.06.026 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A stable P1/P1 finite element for finite strain von Mises elasto-plasticity E. Feulvarcha , J.-C. Rouxa , J.-M. Bergheaua , P. Gillesb a
Univ. Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Etienne cedex 2, France b GEP-INT, 103 rue de Rome, 75017 Paris, France
Abstract The finite element P 1/P 1 is well known for being unsuitable for the simulation of incompressible material flows. In von Mises elasto-plasticity, the volume changes can become negligible when the plastic strain grows higher than the elastic strain. Thus, the material flow is nearly incompressible even if a small volumic elastic strain persists. In this context, the finite element P 1/P 1 leads to pressure oscillations which need to be addressed for achieving satisfactory solutions. The aim of this work is to propose a stabilized formulation without introducing new degrees of freedom or stabilization parameters as for sub-grid scale techniques. Unlike to the standard SUPG method, the stabilization depends on the time step. Examples show that a first order accuracy can be obtained for the pressure and the approach developed can be well suited for cyclic loadings. Keywords: finite element method, mixed formulation, plasticity, stabilization 1. Introduction For von Mises elasto-plastic transformations, the plastic strain can become much higher than the elastic strain and the von Mises criterion imposes a volumic plastic strain equal to zero. Therefore, the material flow can become quasi-incompressible even if a small volumic elastic strain persists. The computation of such a behavior requires treatments similar to the ones used for incompressible material flows. In this context, the finite element method faces a numerical difficulty. The discretization of the kinematics must avoid volumic locking phenomena, while ensuring the uniqueness of the approximation. Locking appears when the kinematic field unable to satisfy both the Email address: [email protected] (E. Feulvarch) Preprint submitted to Elsevier
June 6, 2017
momentum balance and the plastic incompressibility equation. The incompressibility condition can be dealt with by projection methods described in [1] or by penalty techniques which requires numerical recipes such as reduced integration [2] or the B-bar method [3]. Unfortunately, reduced integration cannot be applied to linear tetrahedral elements of type P 1 since the number of integration points is always greater than the number of degrees of freedom associated to the kinematics. However, tetrahedral meshes have the advantage that they can be generated from all kinds of geometries. The quadratic element P 2 may be used for incompressible material flows but it leads to a significant increase of the size of the discrete problem to be solved. The incompressibility can also be taken into account by adding new physical unknowns such as the pressure, as explained in [4]. In this case, the LBB condition (Ladyzhenskaya-Babuˇska-Brezzi) [5, 6] ensures the compatibility of the velocity and the pressure approximations [7] and allows establishing an error estimate. In the literature, several tetrahedral finite elements satisfy this condition such as the elements P 2/P 1 or P 1 + /P 1 [8]. But additional nodes are required for the approximation of the velocity compared to the finite element P 1/P 1 which is not LBB-stable [9]. Indeed, this element is not sensitive to volumetric locking phenomena but it does not ensure the uniqueness of the pressure approximation for perfectly incompressible material flows. However, the recent simulation works of Al Akhrass et al. [10] and Feulvarch et al. [11] show that the same interpolation can be used for the kinematics and the pressure in von Mises elasto-plasticity. Unfortunately, numerical experiments exhibit pressure oscillations which have to be addressed in order to achieve satisfactory solutions. The aim of this work is to propose a new stable finite element P 1/P 1. The first section of this paper is dedicated to the mixed finite element formulation. In the second section, the stabilized formulation is detailed in a Lagrangian approach. It does not require new degrees of freedom or stabilization parameters compared respectively to the OSS sub-grid scale method [12, 13] and the ASGS sub-grid scale method [14]. Unlike to the standard SUPG method [15, 16], the added stabilization term depends on the time step and not only on the spatial discretization. Moreover, the numerical integration is performed using only one centered point. In the last section, examples show that a first order accuracy can be obtained for the pressure approximation in von Mises elasto-plasticity and the approach developed can be well suited for cyclic loadings.
2
2. Finite element modeling 2.1. Strong formulation Within a Lagrangian framework, the momentum in a domain Ω of boundary ∂Ω can be formulated as follows: Find v and p defined on Ω × [0, T ] verifying the initial-boundary value problem defined by ∇ · S + ∇p = 0 in Ω (1a) ∇ · v − 1 p0 = 0 in Ω (1b) α (P ) σ .n = T(p) on ΓT (1c) (p) v=v on Γv (1d)
where p0 = ∂p ; α is the bulk modulus; n is the unit outward normal vector to ∂t the boundary; T(p) is a prescribed “stress vector”; v(p) is a prescribed value of the velocity v; S denotes the deviator part of the Cauchy stress tensor σ. Furthermore, ΓT ∪ Γv = ∂Ω and ΓT ∩ Γv = ∅.
2.2. Weak formulation The weak formulation of the problem (P ) is classically obtained by multiplying Equation (1a) by a weighting function w, Equation (1b) by a weighting function q, and integrating over the domain Ω. Integrating the first equation by parts and accounting for homogeneous Dirichlet boundary conditions for more convenience (see e.g. [17, 18]), one thus obtains the following weak formulation of the problem: Find functions (v, p) ∈ W × Q × [0, T ] such that for all functions (w, q) ∈ W × Q × [0, T ], Z Z Z s T(p) w dS = 0 S : ∇ w dV + p ∇ · w dV − Ω
Ω
Z
Ω
ΓT
q ∇ · v dV −
Z
Ω
1 0 q p dV = 0 α
(2)
with 3 3 W = HΓ1v (Ω) = w ∈ H 1 (Ω) , w = 0 on Γv and Q = L2 (Ω)
∇s w denotes the symmetric part of ∇w.
3
2.3. Spatial integration Following the usual finite element procedure, the finite element P 1/P 1 consists in approximating the velocity v and the pressure p as follows: vh (x) =
N X
vi Ni (x) ; ph (x) =
i=1
N X i=1
h
pi Ni (x) ⇒ p0 (x) =
N X
p0 i Ni (x)
i=1
In these expressions N denotes the number of nodes, vi , pi and p0 i the values of the functions v, p and p0 at node i and Ni (x) the linear shape function associated to this node. Following Galerkin’s standard approach, the functions w and q are taken of the same form. The finite element P 1/P 1 has a very feature for incompressible R interesting 1 h 0h material flows (when α → +∞ i.e. Ω α q p dV → 0). Indeed, it is not sensitive to the volumetric locking phenomenon because the incompressibility equation does not reduce the possible solutions vh to 0: Z
h ∀q , q h ∇ · vh dV = q h [B] {vh } = 0 ⇒ [B] {vh } = {0} 6⇒ {vh } = {0} Ω
where {q h } ≡ (qi )1≤i≤N and {vh } ≡ (vik )1≤i≤N ; 1≤k≤3 denote the vectors of nodal values of the functions q h and vh ; [B] ≡ (Bi,jk )1≤i,j≤N ; 1≤k≤3 is a matrix defined by Z ∂Nj Bi,jk ≡ Ni dV ∂xk Ω
with rank (B) ≤ N < 3 N .
Unfortunatly, the finite element P 1/P 1 is known for being unable to compute incompressible material flows because it does not satisfy the LBB condition given by R h q ∇ · wh dV inf sup Ω h > β0 > 0 q h ∈Qh wh ∈W h kw kW kq h kQ where the constant β0 does not depend on the finite element size h. In pratice, the finite element P 1/P 1 does not guarantee the uniqueness of the pressure approximation ph as explained in [7]. For compressible solid metals, the bulk modulus α does not tend toward +∞ and it has a finite value. In this case, the recent simulation works of Al Akhrass et al. [10] and Feulvarch et al. [11] show that the same interpolation can be used for the kinematics and the pressure in von Mises elasto-plastic transformations. Unfortunately, numerical experiments exhibit pressure oscillations which have 4
to be addressed in order to achieve satisfactory solutions. From the mathematical point of view, the literature does not exhibit any error estimate for the pressure approximation with the finite elementP 1/P 1. The only existing result considers that if the approximation vh , ph ∈ W h × Qh satisfies the LBB condition, the following error estimate holds in linear elasticity [19]: kv − vh kW + kp − ph kQ 6 C inf(qh ∈Qh ,wh ∈W h ) kv − wh kW + kp − q h kQ where the constant C does not depend on the finite element size h.
Fig. 1 presents an example of pressure oscillations obtained with the element P 1/P 1 during the traction of a perfectly plastic axisymmetric specimen. The yield stress σ y and the Poisson’s coefficient ν are respectively taken equal to 300MPa and 0.3. The nodes located at the bottom have a nil vertical displacement while the nodes at the top have a vertical displacement going from 0mm to 0, 4mm in 0, 1s. The simulation has been performed usR ing the computer code Sysweld [20] with a standard finite-strain updated Lagrangian technique dedicated to large transformations in solids [4, 21].
Figure 1: Distribution of pressure (MPa) at the center of each finite element P 1/P 1 during the traction of a perfectly plastic axisymmetric specimen (h=0.52mm and ∆t=1.10−2 s).
5
3. Pressure stabilization To avoid spurious pressure modes, Brezzi andRPitkaranta [22] have first proposed to add a diffusion term of the form βh2 Ω ∇q h .∇ph dV in the discrete form of equation (2)2 . The main drawback of this approach lies in the choice of β which is a non-negative constant. More recently, subgrid scale methods have been developed. One can cite the OSS method [12, 13] or the ASGS method [14] which consist respectively in adding new degrees of freedom or calibrated stabilization parameters. Feulvarch et al. [11] proposed to add an anisotropic diffusion term for the analysis of compressible solid metals (α 9 +∞). This approach seems to give satisfactory results when σ y is nil. Unfortunately, experience shows that the computed stresses are not satisfactory when σ y is not equal to zero. In this work, we propose to compute an anisotropic perturbation term for the elements where the plastic strain varies if they are surrounded by elements which have also a variation of their plastic strain as shown in Fig. 2. In this way, the stabilization is only considered in the quasi-incompressible zone because of the von Mises criterion. Thus, the stabilized problem can be formulated as follows: Find functions vh , ph such that for all functions wh , q h , Z Z Z s h h h S : ∇ w dV + p ∇ · w dV − T(p) wh dS = 0 Ωh
Z
Ωh
q h ∇ · vh dV −
where H = 10 with
Γh T
Ωh
Z
Ωh
1 h 0h q p dV − α
e ΣN i=1
Z
λ(h) =
Z
p
Ωh
∇q h · H · ∇p ∆t dV = 0
10 (D ) Ni dV
Ωh
1 2 (vh )2
(3)
0h
λ(h) vh ⊗ vh
v u 3 uX t (hq vh .xq )2 q=1
hq denotes the length of the considered finite element in the direction xq ; N e is the number of nodes of the element e; Dp is the equivalent plastic strain rate; 10 is an indicator function defined as follows: 1 if x = 0, 10 (x) ≡ 0 if x 6= 0. 6
Figure 2: Selection of the elements to be stabilized versus the equivalent plastic strain rate.
ROne can note that the stabilization term is non-linear and the integration of Ωh α1 q h p0 h dV requires more than one centered integration point for being exactly computed. By considering an implicit updated Lagrangian approach [4, 21], we propose to use a lumping technique [23] and an explicit estimation of H by solving the following discrete form of equation (3)2 at each node i: Z Z Z ∆pi h Ni dV − ∆t ∇Ni .Ht .∇∆ph dV = 0 Ni ∇ · ∆x dV − h h h α Ωt+∆t Ωt+∆t Ωt+∆t where ∆x = vt+∆t ∆t and ∆p = p0 t+∆t ∆t denote respectively the displacement and pressure increments at each computing time t + ∆t. 4. Applications 4.1. Convergence test This application corresponds to the one described in section 2.3. As we have no theoretical error estimate for the pressure approximation, we propose to test the convergence of the stabilized numerical solution ph to the reference solution p plotted in Fig. 3. The field p was built by means of an axisymmetric computation on a 2D mesh made of finite elements of type Q1 employing the B-bar method [3] (h=4.10−2 mm and ∆t=1.10−5 s). Fig. 3 shows that there is no spurious pressure oscillation for h=0.52mm and ∆t=1.10−2 s. In Fig. 4, the order of convergence α is identified by means of the slope of the linear regression on the log-log plot of kp−ph kQ for several time steps (kp − ph kQ ≈ C hα ). One can note that a first order accuracy is obtained regardless of the time step ∆t.
7
s (a)
(b)
Figure 3: Distribution of pressure (MPa) at the center of each finite element: (a) 2D reference solution (h=4.10−2 mm and ∆t=1.10−5 s) - (b) 3D stabilized solution (h=0.52mm and ∆t=1.10−2 s).
Figure 4: Convergence plot of the stabilized numerical solution ph .
8
4.2. Cyclic loading test In this example, we consider a complex elasto-plastic behavior and cyclic loading to show the ability of the stabilized formulation to simulate the complex evolution of a component made of S355 steel. Esnault et al. [24] identified parameters of constitutive equations by considering isotropic and nonlinear kinematic hardening as follows: f = J2 2 (σ p− χ) − Rp χ˙ = 3 c D − γ χ D Rt R = R0 + (Rmax − R0 ) 1 − e−b ( 0 Dp dτ )
where f denotes the yield function; J2 is the second invariant of the deviator; χ˙ is the timeR derivative of the kinematic variable χ; D p is the plastic strain t rate tensor; 0 Dp dτ corresponds to the cumulated equivalent plastic strain; R0 , Rmax , b, c and γ are 5 parameters given in Tab. 1. E (MPa) ν c (MPa) γ R0 (MPa) Rmax (MPa) b
210 000 0.3 46 658 242 346 219 343
Table 1: Coefficients of constitutive equations for S355 steel [24].
Considering such constitutive equations, it is possible to model a softening effect during cyclic loading composed of asymmetric cycles of imposed displacements (see [25, 26]). During softening, the mean stress tends to cyclically relax toward zero. As an example, we consider the cyclic asymmetric linear evolution of the vertical displacement for nodes located at the top of the sample (see Fig. 1). During a cycle, the imposed vertical displacement increases from 0mm to 0,04mm in 0,2s and then, it decreases to 0mm in 0,2s. As shown in Fig. 5 and Fig. 6 after 12,5 cycles, the results are in good agreement with the reference solution built by means of the same 2D mesh (h=4.10−2 mm) and the same time step (∆t=1.10−5 s) as for the previous application. One can note that the results get closer to the reference solution with the refinement of the discretization. As far as the pressure is concerned, no spurious oscillations are observed. Fig. 7 shows the evolution of the vertical stress at the top of the sample. The softening effect clearly appears 9
cycle after cycle and the stationary state is almost reached at 5s for the 3D solution as for the reference solution. 5. Conclusion A stabilized formulation has been proposed for the finite element P1/P1 for the finite-strain analysis of compressible solid metals in the context of von Mises elasto-plastic transformations. The problem to solve is symmetric and the numerical integration is performed using only one centered integration point without adding new degrees of freedom or stabilization parameters to calibrate. Experience shows that a first order accuracy can be obtained for the pressure approximation. Moreover, the results are very satisfactory in the case of cyclic loading. In the future, this work will be extended to the use of the Proper Generalized Decomposition technique for quickly determining the elastoplastic states resulting from cyclic loadings [27].
10
(a) h=0, 52mm-∆t=1.10−2 s (b) h=0, 22mm-∆t=1.10−2 s (c) h=0, 22mm-∆t=1.10−3 s
(d) reference Figure 5: Distribution of pressure (MPa) at the center of each finite element at time t = 5s.
11
(a) h=0, 52mm-∆t=1.10−2 s (b) h=0, 22mm-∆t=1.10−2 s (c) h=0, 22mm-∆t=1.10−3 s
(d) reference Figure 6: Distribution of von Mises stress (MPa) at the center of each finite element at time t = 5s
12
Figure 7: Evolution of the vertical stress at the top of the sample (see Fig. 1).
References [1] Gresho P. M., Sani R. L., Incompressible flow and the finite element method, Wiley, 2000. [2] Reddy J.N., Gartling D.K., The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, 2000. [3] Hughes T.J.R., The Finite Element Method, Dover Publications, Inc. New York, 2000. [4] Bathe K.J., Finite Element Procedure. Prentice Hall, 1996. [5] Babuˇska I., The finite element method with Lagrangien multipliers, Numer Math., vol. 20,pp. 179-192, 1973. [6] Brezzi F., On the existence, uniqueness and approximation of saddlepoint problems arising from Lagrangian multipliers, RAIRO, Anal. Numer, R2, pp. 129-151, 1974. [7] Ern A., Guermond J.-L., Theory and Practice of Finite Elements, Springer Series in Applied Mathematical Sciences, Vol. 159, SpringerVerlag, New York, 2004. [8] Boffi D., Brezzi F., Fortin M., Finite elements for the Stokes problem - Mixed Finite Elements, Compatibility Conditions, and Applications Lecture Notes in Mathematics, vol. 1939, pp. 45-100, 2008. 13
[9] E. Onate, P. Nadukandi, S. R. Idelsohn, P1/P0+ elements for incompressible flows with discontinuous material properties, Comput. Methods Appl. Mech. Engrg. 271 (2014) 185-209. [10] Al Akhrass D., Bruchon J., Drapier S., Fayolle S., Integrating a logarithmic-strain based hyperelastic formulation into a three-field mixed finite element formulation to deal with incompressibility in finite-strain elastoplasticity, Finite Elements in Analysis and Design, Vol. 86, (2014), 61-70. [11] Feulvarch E., Amin El Sayed H., Roux J.-C., Bergheau J.-M., A stabilized P1/P1 finite element for the mechanical analysis of solid metals, Int J Mater Form (2015). doi:10.1007/s12289-015-1252-9. [12] Agelet de Saracibar C., Chiumenti M., Valverde Q., Cervera M., On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity, Comput. Methods Appl. Mech. Engrg., Vol. 195, (2006), 1224-1251. [13] Agelet de Saracibar C. , Chiumenti M., Cervera M., Dialami N., Seret A., Computational Modeling and Sub-Grid Scale Stabilization of Incompressibility and Convection in the Numerical Simulation of Friction Stir Welding Processes, Archives of Computational Methods in Engineering, Vol. 21, (2014), 3-37. [14] Cervera M., Chiumenti M., Benedetti L., Codina R., Mixed stabilized finite element methods in nonlinear solid mechanics. Part III: Compressible and incompressible plasticity, Comput. Methods Appl. Mech. Engrg., Vol. 285, (2015), 752-775. [15] J. Donea, A. Huerta, Finite Element Method for Flow Problems, J. Wiley and Sons, 2003. [16] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu, The Finite Element Method for Fluid Dynamics, Elsevier, 2006. [17] Allaire G., Numerical Analysis and Optimization, Oxford University Press, 2007. [18] Brezis, H., Functional analysis, Sobolev spaces and partial differential equations, Universitext, Springer, 2011. [19] Brezzi, F., Fortin, M., Mixed and hybrid finite element methods, Springer series in computational mathematics (15), Springer Verlag, 1991. 14
[20] ESI Group, Users manual, 2017. [21] Belytschko T., Liu W.K., Moran B., Nonlinear Finite Elements for Continua and Structures, Wiley, 2000. [22] Brezzi F., Pitkaranta J., On the stabilization of finite element approximations of the stokes equations, in Efficient solutions of elliptic systems, Notes on Numerical Fluid Mechanics, Vol. 10 (W. Hackbush, ed.), Braunschweig, Wiesbaden, (1984). [23] E. Feulvarch, J.M. Bergheau, J.B. Leblond, An implicit finite element algorithm for the simulation of diffusion with phase changes in solids, International Journal for Numerical Methods in Engineering, 78 (2009) 1492-1512. [24] J.B. Esnault, V. Doquet, P. Massin, A three-dimensional analysis of fatigue crack paths in thin metallic sheets, International Journal of Fatigue, 62 (2014) 119-132. [25] H. Burlet and G. Cailletaud, Modeling of cyclic plasticity in finite element codes in Desai, C.S. (Ed.), 2nd Int. Conf. On Constitutive Laws for Engineering Materials: Theory and Applications. Elsevier, Tucson, 1157-1164, 1987. [26] J.-L. Chaboche, A review of some plasticity and viscoplasticity constitutive theories, International Journal of Plasticity, 24, 1642-1693, 2008. [27] J.-M. Bergheau, S. Zuchiatti, J.-C. Roux, E. Feulvarch, S. Tissot, G. Perrin, The Proper Generalized Decomposition as a spacetime integrator for elastoplastic problems, C. R. Mecanique 344 (2016) 759768.
15
Highlights : • A stabilized formulation is proposed for finite strain von Mises elastoplasticity. • The stabilization does not require new degrees of freedom or stabilization parameters. • The added perturbation term depends on the time step unlike the SUPG method. • A first order accuracy is obtained for the pressure approximation. • The new formulation is able to model a softening effect due to a non-linear kinematic hardening.
PII: DOI: Reference:
S0045-7825(17)30169-X http://dx.doi.org/10.1016/j.cma.2017.06.026 CMA 11491
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date : 24 January 2017 Revised date : 6 June 2017 Accepted date : 22 June 2017 Please cite this article as: E. Feulvarch, J.-. Roux, J.-. Bergheau, P. Gilles, A stable P1/P1 finite element for finite strain von Mises elasto-plasticity, Comput. Methods Appl. Mech. Engrg. (2017), http://dx.doi.org/10.1016/j.cma.2017.06.026 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A stable P1/P1 finite element for finite strain von Mises elasto-plasticity E. Feulvarcha , J.-C. Rouxa , J.-M. Bergheaua , P. Gillesb a
Univ. Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Etienne cedex 2, France b GEP-INT, 103 rue de Rome, 75017 Paris, France
Abstract The finite element P 1/P 1 is well known for being unsuitable for the simulation of incompressible material flows. In von Mises elasto-plasticity, the volume changes can become negligible when the plastic strain grows higher than the elastic strain. Thus, the material flow is nearly incompressible even if a small volumic elastic strain persists. In this context, the finite element P 1/P 1 leads to pressure oscillations which need to be addressed for achieving satisfactory solutions. The aim of this work is to propose a stabilized formulation without introducing new degrees of freedom or stabilization parameters as for sub-grid scale techniques. Unlike to the standard SUPG method, the stabilization depends on the time step. Examples show that a first order accuracy can be obtained for the pressure and the approach developed can be well suited for cyclic loadings. Keywords: finite element method, mixed formulation, plasticity, stabilization 1. Introduction For von Mises elasto-plastic transformations, the plastic strain can become much higher than the elastic strain and the von Mises criterion imposes a volumic plastic strain equal to zero. Therefore, the material flow can become quasi-incompressible even if a small volumic elastic strain persists. The computation of such a behavior requires treatments similar to the ones used for incompressible material flows. In this context, the finite element method faces a numerical difficulty. The discretization of the kinematics must avoid volumic locking phenomena, while ensuring the uniqueness of the approximation. Locking appears when the kinematic field unable to satisfy both the Email address: [email protected] (E. Feulvarch) Preprint submitted to Elsevier
June 6, 2017
momentum balance and the plastic incompressibility equation. The incompressibility condition can be dealt with by projection methods described in [1] or by penalty techniques which requires numerical recipes such as reduced integration [2] or the B-bar method [3]. Unfortunately, reduced integration cannot be applied to linear tetrahedral elements of type P 1 since the number of integration points is always greater than the number of degrees of freedom associated to the kinematics. However, tetrahedral meshes have the advantage that they can be generated from all kinds of geometries. The quadratic element P 2 may be used for incompressible material flows but it leads to a significant increase of the size of the discrete problem to be solved. The incompressibility can also be taken into account by adding new physical unknowns such as the pressure, as explained in [4]. In this case, the LBB condition (Ladyzhenskaya-Babuˇska-Brezzi) [5, 6] ensures the compatibility of the velocity and the pressure approximations [7] and allows establishing an error estimate. In the literature, several tetrahedral finite elements satisfy this condition such as the elements P 2/P 1 or P 1 + /P 1 [8]. But additional nodes are required for the approximation of the velocity compared to the finite element P 1/P 1 which is not LBB-stable [9]. Indeed, this element is not sensitive to volumetric locking phenomena but it does not ensure the uniqueness of the pressure approximation for perfectly incompressible material flows. However, the recent simulation works of Al Akhrass et al. [10] and Feulvarch et al. [11] show that the same interpolation can be used for the kinematics and the pressure in von Mises elasto-plasticity. Unfortunately, numerical experiments exhibit pressure oscillations which have to be addressed in order to achieve satisfactory solutions. The aim of this work is to propose a new stable finite element P 1/P 1. The first section of this paper is dedicated to the mixed finite element formulation. In the second section, the stabilized formulation is detailed in a Lagrangian approach. It does not require new degrees of freedom or stabilization parameters compared respectively to the OSS sub-grid scale method [12, 13] and the ASGS sub-grid scale method [14]. Unlike to the standard SUPG method [15, 16], the added stabilization term depends on the time step and not only on the spatial discretization. Moreover, the numerical integration is performed using only one centered point. In the last section, examples show that a first order accuracy can be obtained for the pressure approximation in von Mises elasto-plasticity and the approach developed can be well suited for cyclic loadings.
2
2. Finite element modeling 2.1. Strong formulation Within a Lagrangian framework, the momentum in a domain Ω of boundary ∂Ω can be formulated as follows: Find v and p defined on Ω × [0, T ] verifying the initial-boundary value problem defined by ∇ · S + ∇p = 0 in Ω (1a) ∇ · v − 1 p0 = 0 in Ω (1b) α (P ) σ .n = T(p) on ΓT (1c) (p) v=v on Γv (1d)
where p0 = ∂p ; α is the bulk modulus; n is the unit outward normal vector to ∂t the boundary; T(p) is a prescribed “stress vector”; v(p) is a prescribed value of the velocity v; S denotes the deviator part of the Cauchy stress tensor σ. Furthermore, ΓT ∪ Γv = ∂Ω and ΓT ∩ Γv = ∅.
2.2. Weak formulation The weak formulation of the problem (P ) is classically obtained by multiplying Equation (1a) by a weighting function w, Equation (1b) by a weighting function q, and integrating over the domain Ω. Integrating the first equation by parts and accounting for homogeneous Dirichlet boundary conditions for more convenience (see e.g. [17, 18]), one thus obtains the following weak formulation of the problem: Find functions (v, p) ∈ W × Q × [0, T ] such that for all functions (w, q) ∈ W × Q × [0, T ], Z Z Z s T(p) w dS = 0 S : ∇ w dV + p ∇ · w dV − Ω
Ω
Z
Ω
ΓT
q ∇ · v dV −
Z
Ω
1 0 q p dV = 0 α
(2)
with 3 3 W = HΓ1v (Ω) = w ∈ H 1 (Ω) , w = 0 on Γv and Q = L2 (Ω)
∇s w denotes the symmetric part of ∇w.
3
2.3. Spatial integration Following the usual finite element procedure, the finite element P 1/P 1 consists in approximating the velocity v and the pressure p as follows: vh (x) =
N X
vi Ni (x) ; ph (x) =
i=1
N X i=1
h
pi Ni (x) ⇒ p0 (x) =
N X
p0 i Ni (x)
i=1
In these expressions N denotes the number of nodes, vi , pi and p0 i the values of the functions v, p and p0 at node i and Ni (x) the linear shape function associated to this node. Following Galerkin’s standard approach, the functions w and q are taken of the same form. The finite element P 1/P 1 has a very feature for incompressible R interesting 1 h 0h material flows (when α → +∞ i.e. Ω α q p dV → 0). Indeed, it is not sensitive to the volumetric locking phenomenon because the incompressibility equation does not reduce the possible solutions vh to 0: Z
h ∀q , q h ∇ · vh dV = q h [B] {vh } = 0 ⇒ [B] {vh } = {0} 6⇒ {vh } = {0} Ω
where {q h } ≡ (qi )1≤i≤N and {vh } ≡ (vik )1≤i≤N ; 1≤k≤3 denote the vectors of nodal values of the functions q h and vh ; [B] ≡ (Bi,jk )1≤i,j≤N ; 1≤k≤3 is a matrix defined by Z ∂Nj Bi,jk ≡ Ni dV ∂xk Ω
with rank (B) ≤ N < 3 N .
Unfortunatly, the finite element P 1/P 1 is known for being unable to compute incompressible material flows because it does not satisfy the LBB condition given by R h q ∇ · wh dV inf sup Ω h > β0 > 0 q h ∈Qh wh ∈W h kw kW kq h kQ where the constant β0 does not depend on the finite element size h. In pratice, the finite element P 1/P 1 does not guarantee the uniqueness of the pressure approximation ph as explained in [7]. For compressible solid metals, the bulk modulus α does not tend toward +∞ and it has a finite value. In this case, the recent simulation works of Al Akhrass et al. [10] and Feulvarch et al. [11] show that the same interpolation can be used for the kinematics and the pressure in von Mises elasto-plastic transformations. Unfortunately, numerical experiments exhibit pressure oscillations which have 4
to be addressed in order to achieve satisfactory solutions. From the mathematical point of view, the literature does not exhibit any error estimate for the pressure approximation with the finite elementP 1/P 1. The only existing result considers that if the approximation vh , ph ∈ W h × Qh satisfies the LBB condition, the following error estimate holds in linear elasticity [19]: kv − vh kW + kp − ph kQ 6 C inf(qh ∈Qh ,wh ∈W h ) kv − wh kW + kp − q h kQ where the constant C does not depend on the finite element size h.
Fig. 1 presents an example of pressure oscillations obtained with the element P 1/P 1 during the traction of a perfectly plastic axisymmetric specimen. The yield stress σ y and the Poisson’s coefficient ν are respectively taken equal to 300MPa and 0.3. The nodes located at the bottom have a nil vertical displacement while the nodes at the top have a vertical displacement going from 0mm to 0, 4mm in 0, 1s. The simulation has been performed usR ing the computer code Sysweld [20] with a standard finite-strain updated Lagrangian technique dedicated to large transformations in solids [4, 21].
Figure 1: Distribution of pressure (MPa) at the center of each finite element P 1/P 1 during the traction of a perfectly plastic axisymmetric specimen (h=0.52mm and ∆t=1.10−2 s).
5
3. Pressure stabilization To avoid spurious pressure modes, Brezzi andRPitkaranta [22] have first proposed to add a diffusion term of the form βh2 Ω ∇q h .∇ph dV in the discrete form of equation (2)2 . The main drawback of this approach lies in the choice of β which is a non-negative constant. More recently, subgrid scale methods have been developed. One can cite the OSS method [12, 13] or the ASGS method [14] which consist respectively in adding new degrees of freedom or calibrated stabilization parameters. Feulvarch et al. [11] proposed to add an anisotropic diffusion term for the analysis of compressible solid metals (α 9 +∞). This approach seems to give satisfactory results when σ y is nil. Unfortunately, experience shows that the computed stresses are not satisfactory when σ y is not equal to zero. In this work, we propose to compute an anisotropic perturbation term for the elements where the plastic strain varies if they are surrounded by elements which have also a variation of their plastic strain as shown in Fig. 2. In this way, the stabilization is only considered in the quasi-incompressible zone because of the von Mises criterion. Thus, the stabilized problem can be formulated as follows: Find functions vh , ph such that for all functions wh , q h , Z Z Z s h h h S : ∇ w dV + p ∇ · w dV − T(p) wh dS = 0 Ωh
Z
Ωh
q h ∇ · vh dV −
where H = 10 with
Γh T
Ωh
Z
Ωh
1 h 0h q p dV − α
e ΣN i=1
Z
λ(h) =
Z
p
Ωh
∇q h · H · ∇p ∆t dV = 0
10 (D ) Ni dV
Ωh
1 2 (vh )2
(3)
0h
λ(h) vh ⊗ vh
v u 3 uX t (hq vh .xq )2 q=1
hq denotes the length of the considered finite element in the direction xq ; N e is the number of nodes of the element e; Dp is the equivalent plastic strain rate; 10 is an indicator function defined as follows: 1 if x = 0, 10 (x) ≡ 0 if x 6= 0. 6
Figure 2: Selection of the elements to be stabilized versus the equivalent plastic strain rate.
ROne can note that the stabilization term is non-linear and the integration of Ωh α1 q h p0 h dV requires more than one centered integration point for being exactly computed. By considering an implicit updated Lagrangian approach [4, 21], we propose to use a lumping technique [23] and an explicit estimation of H by solving the following discrete form of equation (3)2 at each node i: Z Z Z ∆pi h Ni dV − ∆t ∇Ni .Ht .∇∆ph dV = 0 Ni ∇ · ∆x dV − h h h α Ωt+∆t Ωt+∆t Ωt+∆t where ∆x = vt+∆t ∆t and ∆p = p0 t+∆t ∆t denote respectively the displacement and pressure increments at each computing time t + ∆t. 4. Applications 4.1. Convergence test This application corresponds to the one described in section 2.3. As we have no theoretical error estimate for the pressure approximation, we propose to test the convergence of the stabilized numerical solution ph to the reference solution p plotted in Fig. 3. The field p was built by means of an axisymmetric computation on a 2D mesh made of finite elements of type Q1 employing the B-bar method [3] (h=4.10−2 mm and ∆t=1.10−5 s). Fig. 3 shows that there is no spurious pressure oscillation for h=0.52mm and ∆t=1.10−2 s. In Fig. 4, the order of convergence α is identified by means of the slope of the linear regression on the log-log plot of kp−ph kQ for several time steps (kp − ph kQ ≈ C hα ). One can note that a first order accuracy is obtained regardless of the time step ∆t.
7
s (a)
(b)
Figure 3: Distribution of pressure (MPa) at the center of each finite element: (a) 2D reference solution (h=4.10−2 mm and ∆t=1.10−5 s) - (b) 3D stabilized solution (h=0.52mm and ∆t=1.10−2 s).
Figure 4: Convergence plot of the stabilized numerical solution ph .
8
4.2. Cyclic loading test In this example, we consider a complex elasto-plastic behavior and cyclic loading to show the ability of the stabilized formulation to simulate the complex evolution of a component made of S355 steel. Esnault et al. [24] identified parameters of constitutive equations by considering isotropic and nonlinear kinematic hardening as follows: f = J2 2 (σ p− χ) − Rp χ˙ = 3 c D − γ χ D Rt R = R0 + (Rmax − R0 ) 1 − e−b ( 0 Dp dτ )
where f denotes the yield function; J2 is the second invariant of the deviator; χ˙ is the timeR derivative of the kinematic variable χ; D p is the plastic strain t rate tensor; 0 Dp dτ corresponds to the cumulated equivalent plastic strain; R0 , Rmax , b, c and γ are 5 parameters given in Tab. 1. E (MPa) ν c (MPa) γ R0 (MPa) Rmax (MPa) b
210 000 0.3 46 658 242 346 219 343
Table 1: Coefficients of constitutive equations for S355 steel [24].
Considering such constitutive equations, it is possible to model a softening effect during cyclic loading composed of asymmetric cycles of imposed displacements (see [25, 26]). During softening, the mean stress tends to cyclically relax toward zero. As an example, we consider the cyclic asymmetric linear evolution of the vertical displacement for nodes located at the top of the sample (see Fig. 1). During a cycle, the imposed vertical displacement increases from 0mm to 0,04mm in 0,2s and then, it decreases to 0mm in 0,2s. As shown in Fig. 5 and Fig. 6 after 12,5 cycles, the results are in good agreement with the reference solution built by means of the same 2D mesh (h=4.10−2 mm) and the same time step (∆t=1.10−5 s) as for the previous application. One can note that the results get closer to the reference solution with the refinement of the discretization. As far as the pressure is concerned, no spurious oscillations are observed. Fig. 7 shows the evolution of the vertical stress at the top of the sample. The softening effect clearly appears 9
cycle after cycle and the stationary state is almost reached at 5s for the 3D solution as for the reference solution. 5. Conclusion A stabilized formulation has been proposed for the finite element P1/P1 for the finite-strain analysis of compressible solid metals in the context of von Mises elasto-plastic transformations. The problem to solve is symmetric and the numerical integration is performed using only one centered integration point without adding new degrees of freedom or stabilization parameters to calibrate. Experience shows that a first order accuracy can be obtained for the pressure approximation. Moreover, the results are very satisfactory in the case of cyclic loading. In the future, this work will be extended to the use of the Proper Generalized Decomposition technique for quickly determining the elastoplastic states resulting from cyclic loadings [27].
10
(a) h=0, 52mm-∆t=1.10−2 s (b) h=0, 22mm-∆t=1.10−2 s (c) h=0, 22mm-∆t=1.10−3 s
(d) reference Figure 5: Distribution of pressure (MPa) at the center of each finite element at time t = 5s.
11
(a) h=0, 52mm-∆t=1.10−2 s (b) h=0, 22mm-∆t=1.10−2 s (c) h=0, 22mm-∆t=1.10−3 s
(d) reference Figure 6: Distribution of von Mises stress (MPa) at the center of each finite element at time t = 5s
12
Figure 7: Evolution of the vertical stress at the top of the sample (see Fig. 1).
References [1] Gresho P. M., Sani R. L., Incompressible flow and the finite element method, Wiley, 2000. [2] Reddy J.N., Gartling D.K., The Finite Element Method in Heat Transfer and Fluid Dynamics, CRC Press, 2000. [3] Hughes T.J.R., The Finite Element Method, Dover Publications, Inc. New York, 2000. [4] Bathe K.J., Finite Element Procedure. Prentice Hall, 1996. [5] Babuˇska I., The finite element method with Lagrangien multipliers, Numer Math., vol. 20,pp. 179-192, 1973. [6] Brezzi F., On the existence, uniqueness and approximation of saddlepoint problems arising from Lagrangian multipliers, RAIRO, Anal. Numer, R2, pp. 129-151, 1974. [7] Ern A., Guermond J.-L., Theory and Practice of Finite Elements, Springer Series in Applied Mathematical Sciences, Vol. 159, SpringerVerlag, New York, 2004. [8] Boffi D., Brezzi F., Fortin M., Finite elements for the Stokes problem - Mixed Finite Elements, Compatibility Conditions, and Applications Lecture Notes in Mathematics, vol. 1939, pp. 45-100, 2008. 13
[9] E. Onate, P. Nadukandi, S. R. Idelsohn, P1/P0+ elements for incompressible flows with discontinuous material properties, Comput. Methods Appl. Mech. Engrg. 271 (2014) 185-209. [10] Al Akhrass D., Bruchon J., Drapier S., Fayolle S., Integrating a logarithmic-strain based hyperelastic formulation into a three-field mixed finite element formulation to deal with incompressibility in finite-strain elastoplasticity, Finite Elements in Analysis and Design, Vol. 86, (2014), 61-70. [11] Feulvarch E., Amin El Sayed H., Roux J.-C., Bergheau J.-M., A stabilized P1/P1 finite element for the mechanical analysis of solid metals, Int J Mater Form (2015). doi:10.1007/s12289-015-1252-9. [12] Agelet de Saracibar C., Chiumenti M., Valverde Q., Cervera M., On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity, Comput. Methods Appl. Mech. Engrg., Vol. 195, (2006), 1224-1251. [13] Agelet de Saracibar C. , Chiumenti M., Cervera M., Dialami N., Seret A., Computational Modeling and Sub-Grid Scale Stabilization of Incompressibility and Convection in the Numerical Simulation of Friction Stir Welding Processes, Archives of Computational Methods in Engineering, Vol. 21, (2014), 3-37. [14] Cervera M., Chiumenti M., Benedetti L., Codina R., Mixed stabilized finite element methods in nonlinear solid mechanics. Part III: Compressible and incompressible plasticity, Comput. Methods Appl. Mech. Engrg., Vol. 285, (2015), 752-775. [15] J. Donea, A. Huerta, Finite Element Method for Flow Problems, J. Wiley and Sons, 2003. [16] O.C. Zienkiewicz, R.L. Taylor, P. Nithiarasu, The Finite Element Method for Fluid Dynamics, Elsevier, 2006. [17] Allaire G., Numerical Analysis and Optimization, Oxford University Press, 2007. [18] Brezis, H., Functional analysis, Sobolev spaces and partial differential equations, Universitext, Springer, 2011. [19] Brezzi, F., Fortin, M., Mixed and hybrid finite element methods, Springer series in computational mathematics (15), Springer Verlag, 1991. 14
[20] ESI Group, Users manual, 2017. [21] Belytschko T., Liu W.K., Moran B., Nonlinear Finite Elements for Continua and Structures, Wiley, 2000. [22] Brezzi F., Pitkaranta J., On the stabilization of finite element approximations of the stokes equations, in Efficient solutions of elliptic systems, Notes on Numerical Fluid Mechanics, Vol. 10 (W. Hackbush, ed.), Braunschweig, Wiesbaden, (1984). [23] E. Feulvarch, J.M. Bergheau, J.B. Leblond, An implicit finite element algorithm for the simulation of diffusion with phase changes in solids, International Journal for Numerical Methods in Engineering, 78 (2009) 1492-1512. [24] J.B. Esnault, V. Doquet, P. Massin, A three-dimensional analysis of fatigue crack paths in thin metallic sheets, International Journal of Fatigue, 62 (2014) 119-132. [25] H. Burlet and G. Cailletaud, Modeling of cyclic plasticity in finite element codes in Desai, C.S. (Ed.), 2nd Int. Conf. On Constitutive Laws for Engineering Materials: Theory and Applications. Elsevier, Tucson, 1157-1164, 1987. [26] J.-L. Chaboche, A review of some plasticity and viscoplasticity constitutive theories, International Journal of Plasticity, 24, 1642-1693, 2008. [27] J.-M. Bergheau, S. Zuchiatti, J.-C. Roux, E. Feulvarch, S. Tissot, G. Perrin, The Proper Generalized Decomposition as a spacetime integrator for elastoplastic problems, C. R. Mecanique 344 (2016) 759768.
15
Highlights : • A stabilized formulation is proposed for finite strain von Mises elastoplasticity. • The stabilization does not require new degrees of freedom or stabilization parameters. • The added perturbation term depends on the time step unlike the SUPG method. • A first order accuracy is obtained for the pressure approximation. • The new formulation is able to model a softening effect due to a non-linear kinematic hardening.