An axisymmetric nodal averaged finite element

Abstract

A nodal averaging technique which was earlier used for plane strain and three-dimensional problems is extended to include the axisymmetric one. Based on the virtual work principle, an expression for nodal force is found. In turn, a nodal force variation yields a stiffness matrix that proves to be non-symmetrical. But, cumbersome non-symmetrical terms can be rejected without the loss of Newton-Raphson iterations convergence. An approximate formula of volume for a ring of triangular profile is exploited in order to simplify program codes and also to accelerate calculations. The proposed finite element is intended primarily for quasistatic problems and large irreversible strain i.e. for metal forming analysis. As a test problem, deep rolling of a steel rod is studied.

Keywords
finite element; nodal averaging; axisymmetric problem; large strain; metal forming

1 INTRODUCTION

The finite element method (FEM) is widely used in science and engineering for numerical analysis of problems which beggar analytical methods. The specific character of solid mechanics tasks is that (except for porous media) a body shape can be varied greatly while body’s volume alters insignificantly. But, the approximating functions proposed by classical FEM (see for example the textbook Zienkiewicz (1971) ) do not comply with the incompressibility requirement. At first sight, this does not induce any drawback because the convergence theorems have been proven for many FEM calculation schemes. Nevertheless, in practice (mainly when a personal computer is used) an approximate solution can be found to be too far from an exact one and such tasks were really discussed in last decades. This circumstance forced to revise all the classical approximations of tensor fields. Many new approaches have been developed but they require a special review. Plane strain, three-dimensional, and shell problems have been involved but the axisymmetric one remains out of consideration so far and the present work aim is elimination of this gap.

The first attempt to overcome the volumetric locking difficulties was the application of selective reduced integration schemes to displacement-based nonlinear finite elements (bi-linear, tri-linear, second order, and third order). The deviatoric items of a variational equation are integrated by larger number of Gauss points P₁,…, P_n while the volumetric items – by less number of points Q₁,…,Q_m ( m<n). In general, it is impossible to meet the incompressibility requirements at all the points P_i but it is possible to do it at points Q_i. This method suffers from many shortcomings: sizable computer memory is needed to store the state parameters such as stress, strain, damage, etc. at Gauss points; problems appear when transferring data from one finite element mesh to another during the remeshing procedure; and so on. Nevertheless, the reduced-integration technique came into further development in Reese (2002 , 2005 ) but is not popular now.

Mixed formulations (see for example Simo et al. (1985) ) are numerous and widespread. They are employed in many commercial software products: ANSYS, NASTRAN, ABAQUS, SUPERFORGE, etc. Degrees of freedom comprise both nodal coordinates and nodal stresses. As a variant, only the hydrostatic component of the stress tensor is used in Sussman and Bathe (1987) . A finite element matrix appears to be sign-indefinite and may be ill-conditioned such that some preconditioner may be required. The convergence analysis of mixed formulations is presented in monograph by Brenner and Scott (2008) . The comparative analysis of a mixed constant pressure hexahedral element Q1P0 and a stabilized nodally integrated tetrahedral by Puso and Solberg (2006) had shown that the results were identical in some tests while the benchmark “Cook membrane” demonstrated the advantage of element by Puso and Solberg. The Q1P0 element converges slowly in this problem due to poor bending performance.

The mixed methods came into further development. A second order stabilized Petrov-Galerkin finite element framework was introduced by Lee et al. (2014) . This formulation, written as a system of conservation laws for the momentum and the deformation gradient, yields equal order of convergence for displacements and stresses. In Gil et al. (2014) , the formulation is enhanced for nearly and truly incompressible materials with some novelties. Unfortunately, that work is restricted by hyper-elastic scenarios.

Enhanced-assumed strain approach has produced many treatments – see for example Simo and Armero (1992) , Simo et al. (1993) , Puso (2000) , and Areias et al. (2003) . An attractive variant was presented by Broccardo et al. (2009) where the assumed strain technique is combined with the nodal average strain formulation but only elastic and hyper-elastic scenarios were under consideration.

Nodal integration, or nodal averaging, yields a particle-type method where stress and material history are located exclusively at the nodes and can be employed when using meshless or finite element shape functions. This particle feature is desirable for large deformation settings because it avoids the remapping or advection of the state parameters required in other methods. The nodal integration relies on fewer stress point evaluation than most other methods ( Puso et al. (2008) ). Three classical nodal integration schemes are known: nodal strain method by Beissel and Belytschko (1996) , stabilized conforming nodal integration by Chen et al. (2001) , and nodal averaging by Dohrmann et al. (2000) . The enhanced versions of these schemes are given in Puso et al. (2008) where the stability of nodal integration for both meshless and finite element approximations was investigated.

Meshfree, or meshless, methods do not employ a finite element base functions therefore these methods do not suffer from the mesh distortion. More complex (in comparison to finite element) meshfree base functions provide with smooth solutions but computational cost may appear to be too high. A comparative study of the three meshless approaches was carried out in Quak et al. (2011) and a comprehensive review is given by Cueto and Chinesta (2015) .

Taking into account all the aforecited, one can conclude that there exist many ways to design an axial symmetric finite element involving super-large strain and high gradients of tensor fields. Therefore, to make an appropriate choice, some special requirements are needed. In this work, the following requirements are determinative:

1. 1

   the volumetric locking is damped;

2. 2

   the locking due to bending is damped;

3. 3

   spurious modes are absent;

4. 4

   a finite element passes the “patch test”;

5. 5

   a theoretical framework is simple;

6. 6

   a program code is short.

A finite element by Puso and Solberg (2006) meet all the requirements and its formulation (based on the nodal integration technique – see above) will be extended on the axial symmetric event here. Unfortunately, two shortcomings are inherent to the formulation. First, mesh distortion causes bad performance and a remeshing procedure together with a data transfer procedure (from an old mesh to a new one) are needed is such case. This problem is studied and resolved in Zhang and Dolbow (2017) specially for the nodal averaged stabilized element by Puso and Solberg (2006) . Second, despite exhibiting very good behavior in terms of displacements the formulation tends to exhibit non-physical hydrostatic pressure fluctuations – see Pires et al. (2004) and Gee et al. (2009) . In order to overcome this drawback, a smoothing procedure is applied at the final stage of numerical analysis – see Puso and Solberg (2006) . Smoothed pressures are evaluated by averaging nodal pressures at element centers and then re-averaging these element pressures back out to the nodes. If it is necessary, this procedure is repeated several times in order to smooth the isoline pattern. So, both the shortcomings of the nodal averaged stabilized element can be rectified.

An approximation based on the simplest finite element (the linear triangle) and the nodal averaging procedure over triangles of a nodal complex is proposed in this work. Due to many reasons, finite elements of low order prove to be preferable for non-linear problems and large deformation. Besides, mesh generators provides triangles and tetrahedrons as a rule; more complex figures especially for irregular meshes cause serious problems. The procedure of nodal averaging was first formulated by Dohrmann et al. (2000) and further analyzed by Bonet et al. (2001) . As it turns out, nodal averaging provides weakened constraints such that coarse mesh accuracy can often be achieved for near incompressible materials. It was pointed out also that the formulation was prone to spurious low-energy modes. In order to suppress them, Puso and Solberg (2006) proposed to enter stabilizing terms into a discrete form of variational equation. A number of different formulations have been developed to rectify poor performance of low-order elements under incompressible and near incompressible conditions but the approach by Dohrmann et al. (2000) alleviates both the volumetric locking and locking due to bending. However, some modern methods give the same result as in Puso and Solberg (2006) – see for example Puso et al. (2008) . But, stabilized nodal averaging is the simplest among them from the point of view of both a theoretical framework and a program code. The present work extends that formulation on the axisymmetric problem.

This work has one more aspect associated with an impotent industrial application – numerical simulation of the deep rolling working. Such the treatment creates a near-surface hardened layer and compressive residual stress preventing from crack generation in metallic parts. Traditionally, this simulation is performed in three dimensional setting. But, in some events it is possible to simplify a model to an axisymmetric one and to reduce the computational cost drastically – see Section 3.

2 FINITE ELEMENT FORMULATION

Let an arbitrary triangulation (i.e. N nodes $A_1\mathrm{,...,}{A}_N$ and M triangles $T_1\mathrm{,...,}{T}_M$ ) of a plane region is given. Let a nodal complex $C_n$ be a set of triangles sharing node $A_n$ . Denote by $<n>$ a set of triangle numbers in complex $C_n$ ; by $<np>$ – a set of triangle numbers in intersection of complexes $C_n$ and $C_p$ (i.e. $<np>=<n>\cap <p>$ ); by $[n]$ – a set of vertex numbers in complex $C_n$ ; by $]m[$ – a set of node numbers of triangular $T_m$ . Figure 1 illustrates these denotations.

Figure 1
Nodal complexes.

The idea of stabilized nodal averaging can be demonstrated by the simplest example – plane strain problem with small elastic strain. Let $u=(u^1,u^2)$ be a FEM-approximation of the displacement field. Define average deformation ${\overline{\epsilon }}_n$ at node $A_n$ and average area ${\overline{S}}_n$ as

${\overline{\epsilon }}_n={\displaystyle \sum _{i\in <n>}{w}_i{\epsilon }_i}$ and ${\overline{S}}_n=S_n/3$ (1)

where ${\epsilon }_i$ – deformation on triangle $T_i$ ( ${\epsilon }_i\equiv const$ because $T_i$ is a linear finite element); $w_i=s_i/S_n$ – weight coefficient; $s_i$ – area of $T_i$ ; $S_n$ – area of $C_n$ . Write the discretized variational equation taking into account (1):

${\displaystyle \sum _{n=1}^N{\overline{S}}_n}\delta {\overline{\epsilon }}_n\cdot \cdot (C-\stackrel{⌢}{C})\cdot \cdot {\overline{\epsilon }}_n+{\displaystyle \sum _{m=1}^M{s}_m}\delta {\epsilon }_m\cdot \cdot \stackrel{⌢}{C}\cdot \cdot {\overline{\epsilon }}_m=f\cdot \delta u$ (2)

where the first sum corresponds with nodal complexes while the second one – with triangles; f is an external load; $C$ is the isotropic elasticity tensor and $\stackrel{⌢}{C}$ is the stabilization tensor (on optimal choice of $\stackrel{⌢}{C}$ see Puso and Solberg (2006) ).

It is clear that the influence of the stabilization tensor drops as a finite element mesh is refined but for a coarse mesh tensor $\stackrel{⌢}{C}$ becomes to be necessary. Formally, one can consider that the tessellation procedure follows the meshing one and new finite elements appear while each of them is a one-third of a nodal complex – see Figure 2 where one-thirds of triangles (each of triangles is dissected into three parts of equal area) form a new (nodal averaged) finite element. Note that all the nodes except for the central one are located outside the new finite element.

Figure 2
Nodal finite elements with 7 and 9 nodes.

As for the axisymmetric event, a volume formula for a ring of triangular profile is not so trivial as for triangular or tetrahedron. This results in too cumbersome expressions for both nodal force and stiffness matrix. In the present work, an approximate formula of volume for such a ring is exploited in order to simplify a program code and to accelerate calculations. Turn to the axisymmetric formulation itself. In cylindrical coordinates $y^1=r$ , $y^2=h$ , $y^3=\varphi$ , the Cauchy stress tensor and the small strain one take the form (here, the small displacement field $u=(u_r,u_h)$ at any time interval $[t,t+\Delta t]$ is meant)

$\sigma =\left[\begin{array}{ccc}{\sigma }_r& {\sigma }_{rh}& 0\\ {\sigma }_{rh}& {\sigma }_h& 0\\ 0& 0& {\sigma }_{\varphi }\end{array}\right]$ , $\epsilon =\frac{1}{2}\left[\begin{array}{ccc}2u_{r,r}& u_{r,h}+u_{h,r}& 0\\ u_{r,h}+u_{h,r}& 2u_{h,h}& 0\\ 0& 0& 2u_r/r\end{array}\right]$ (3)

It is necessary to rewrite the variational principle (2) taking into account large irreversible (plastic) strain and the cpecific form (3) of tensors σ and ε (the Euler approach is adopted throughout the work i.e. all tensor fields refer to the current configuration rather than the initial one):

${\displaystyle \sum _{n=1}^N{\overline{V}}_n{\overline{\sigma }}_n\cdot \cdot }\delta {\epsilon }_n+{\displaystyle \sum _{m=1}^M{v}_m{\stackrel{⌢}{\sigma }}_m\cdot \cdot }\delta {\epsilon }_m=f\cdot \delta u$ (4)

where ${\overline{\sigma }}_n={\displaystyle \sum _{i\in <n>}{w}_i{\sigma }_i}$ – average stress at node $A_n$ ;

${\stackrel{⌢}{\sigma }}_m$ – stabilizing stress on triangle $T_m$ ;

${\overline{\epsilon }}_n={\displaystyle \sum _{i\in <n>}{w}_i{\epsilon }_i}$ – average strain at node $A_n$ ;

$w_i=v_i/V_n$ – weight coefficient of $T_i$ in nodal complex $C_n$ ;

${\overline{V}}_n=V_n/3$ – averaged nodal volume;

$V_n={\displaystyle \sum _{m\in <n>}{v}_m}$ – volume of nodal complex $C_n$ ;

$v_m=2\pi r_m{s}_m$ – approximate volume for a ring of triangular profile $T_m$ ;

$r_m=(r_1+r_2+r_3)/3$ – average radius of a ring of triangular profile $T_m$ .

Equation (4) has to be supplemented by a constitutive equation. As far as metals and alloys refer to near rigid-plastic materials, it is convenient to write that equation by means of any corotational rate “ r” (for example, Yaumann’s rate “j ”):

${}^r\sigma =(C-\stackrel{⌢}{C})\cdot \cdot d$ , ${}^r\stackrel{⌢}{\sigma }=\stackrel{⌢}{C}\cdot \cdot d$ (5)

where d is the deformation rate tensor. Therefore, two stress tensors take part in step-by-step numerical analysis: real stress $\sigma$ and stabilizing stress $\stackrel{⌢}{\sigma }$ .

Introduce a FEM-approximation for the displacement field u, its gradient , and circular deformation ${\epsilon }_{33}$:

$u^i=u_n^i{N}_n(y^1,y^2)$ ; $u_{,j}^i=u_n^i{N}_{n,j}(y^1,y^2)(1\le i,j\le 2)$ ; ${\epsilon }_{33}=u{}_r/r=u_n^1{N}_n(y^1,y^2)/y^1$ (6)

where the summation over repeated indexes is implied ; the differentiation sign is separated by a comma; $N_n(y^1,y^2)$ – shape function corresponding to node $A_n$ ; $(u_n^1,u_n^2)$ – nodal displacement of $A_n$ . In the sequel, a shape function will be numerated by a pair of indexes: $N_n^m(y^1,y^2)$ where n corresponds to node $A_n$ and m – to triangle $T_m$ . Define an average nodal gradient as

${\overline{u}}_{,j}^{(n)i}={\displaystyle \sum _{m\in <n>}{w}_m^n{u}_{,j}^{(m)i}(1\le i,j\le 2)}$ (7)

where ${\overline{u}}_{,j}^{(n)i}$ and $u_{,j}^{(m)i}$ correspond to node A_n and triangle T_m respectively. FEM-approximation of (7), based on approximation (6), yields

${\overline{u}}_{,j}^{(n)i}={\displaystyle \sum _{m\in <n>}{w}_m^n{\displaystyle \sum _{p\in ]m[}{N}_{p,j}^m}{u}_p^i={\displaystyle \sum _{p\in [n]}{u}_p^i{\displaystyle \sum _{m\in <pn>}{w}_m^n{N}_{p,j}^m}={\displaystyle \sum _{p\in [n]}{u}_p^i{\overline{B}}_{pj}^n}}(1\le i,j\le 2)}$ (8)

while the summation in the last but one sum runs over the intersection of nodal complexes $C_p$ and $C_n$ (this intersection can be empty, consists of one or two triangles or coincides with the entire nodal complex if $n=p$ ); matrix ${\overline{B}}^{(n)}$ is of the form

${\overline{B}}_{pj}^n={\displaystyle \sum _{m\in <pn>}{w}_m^n{N}_{p,j}^m}$ (9)

Circular deformation ${\epsilon }_{33}$ is calculated according to

${\epsilon }_{33}={\displaystyle \sum _{p\in [n]}{u}_p^1{\tilde{B}}_p^n}$ (10)

where vector ${\tilde{B}}^{(n)}$ is of the form

${\tilde{B}}_p^n={\displaystyle \sum _{m\in <pn>}{w}_m^n{N}_p^m(y^1,y^2)/y^1}$ (11)

Due to the approximate formula for volume $v_m$ , expression (9) yields

${\overline{B}}_{pj}^n={\displaystyle \sum _{m\in <pn>}{w}_m^n{N}_{p,j}^m}={\displaystyle \sum _{m\in <pn>}\frac{v_m}{V_n}}\frac{b_{pj}^m}{D_m}=\frac{2\pi }{V_n}{\displaystyle \sum _{m\in <pn>}{r}_m{s}_m}\frac{b_{pj}^m}{D_m}=\frac{\pi }{V_n}{\overline{b}}_{pj}^n(1\le i,j\le 2)$ (12)

where $D_m=2s_m$ ; ${\overline{b}}_{pj}^n={\displaystyle \sum _{m\in <pn>}{r}_m}{b}_{pj}^m$ ; $b^{(m)}$ is the well-known in FEM matrix of coordinate differences of triangle vertices (three nodes are numerated by 1,2,3 counter clockwise and index m is omitted):

$b=\left[\begin{array}{cc}{y}_2-y_3& x_3-x_2\\ y_3-y_1& x_1-x_3\\ y_1-y_2& x_2-x_1\end{array}\right]$ . (13)

As for expression (11) , taking into account the one-point integration scheme for the classical linear triangle element, function $N_p^m(y^1,y^2)/y^1$ should be evaluated at the center of triangle $T_m$ where this function equals to $1/(3r_m)$ such that

${\tilde{B}}_p^n=\frac{1}{3}{\displaystyle \sum _{m\in <pn>}\frac{v_m}{V_n{r}_m}}=\frac{2\pi }{3V_n}{\displaystyle \sum _{m\in <pn>}{s}_m}=\frac{2\pi }{3V_n}{\tilde{b}}_p^n$ (14)

where ${\tilde{b}}_p^n={\tilde{b}}_n^p={\displaystyle \sum _{m\in <pn>}{s}_m}$ .

Now, return to variational principle (4) . Substituting the displacement gradient approximation (8) and the approximation of circular deformation (10) , find a nodal force at node $A_p$:

$f_p^i={\displaystyle \sum _{n\in [p]}{\overline{V}}_n{\overline{B}}_{pj}^n{\overline{\sigma }}_n^{ij}+{\delta }^{i1}{\displaystyle \sum _{n\in [p]}{\overline{V}}_n{\tilde{B}}_p^n{\overline{\sigma }}_n^{33}+}}{\displaystyle \sum _{m\in <p>}{v}_m{N}_{p,j}^m{\stackrel{⌢}{\sigma }}_m^{ij}+{\delta }^{i1}{\displaystyle \sum _{m\in <p>}{v}_m{\stackrel{⌢}{\sigma }}_m^{33}{N}_p^m/y^1}}(1\le i,j\le 2)$ (15)

where ${\delta }^{ij}$ is Kronecker’s delta-symbol. Substituting (12) and (14) into (15) and taking into account ${\overline{V}}_n=V_n/3$ , obtain

$f_p^i=\frac{\pi }{3}{\displaystyle \sum _{n\in [p]}{\overline{b}}_{pj}^n{\overline{\sigma }}_n^{ij}+{\delta }^{i1}\frac{2\pi }{9}{\displaystyle \sum _{n\in [p]}{\tilde{b}}_p^n{\overline{\sigma }}_n^{33}+}}\pi {\displaystyle \sum _{m\in <p>}{r}_m{b}_{pj}^m{\stackrel{⌢}{\sigma }}_m^{ij}+{\delta }^{i1}\frac{2\pi }{3}{\displaystyle \sum _{m\in <p>}{s}_m{\stackrel{⌢}{\sigma }}_m^{33}}}(1\le i,j\le 2)$ (16)

Point out that this expression is simpler than (15) because all volumes (both nodal and of rings of triangular profile) disappeared due to the successful choice of the approximate formula for volume $v_m$ (see the elucidations to equation (4) ). Finally, a stiffness matrix is obtained by means of the nodal force variation. To this goal, find the stress tensor variation from the constitutive equation (5) with Yaumann’s rate:

$\dot{\sigma }={}^j\sigma -\sigma \cdot \omega +\omega \cdot \sigma =(C-\stackrel{⌢}{C})\cdot \cdot d-\sigma \cdot \omega +\omega \cdot \sigma$ ; (17)

$\dot{\stackrel{⌢}{\sigma }}={}^j\stackrel{⌢}{\sigma }-\stackrel{⌢}{\sigma }\cdot \omega +\omega \cdot \stackrel{⌢}{\sigma }=\stackrel{⌢}{C}\cdot \cdot d-\stackrel{⌢}{\sigma }\cdot \omega +\omega \cdot \stackrel{⌢}{\sigma }$ (18)

where $\omega =(\nabla v-\nabla v^T)/2$ is the vortex tensor. According to (17) , the nodal stress variation at node $A_n$ is

$\delta {\overline{\sigma }}_n^{ij}=(C_{\overline{n}}^{ijkl}-{\stackrel{⌢}{C}}_{\overline{n}}^{ijkl})\delta {\overline{\epsilon }}_{kl}^{\overline{n}}+\frac{1}{2}{\displaystyle \sum _{q\in [n]}\delta u_q^k(-{\sigma }_{\overline{n}}^{kj}{\overline{B}}_{qi}^{\overline{n}}-{\sigma }_{\overline{n}}^{ik}{\overline{B}}_{qj}^{\overline{n}}+{\delta }_k^i{\sigma }_{\overline{n}}^{lj}{\overline{B}}_{ql}^{\overline{n}}+{\delta }_k^j{\sigma }_{\overline{n}}^{li}{\overline{B}}_{ql}^{\overline{n}})}(1\le i,j,k,l\le 2)$ (19)

where $\overline{n}$ denotes “no summation” over index n. Besides, matricies ${\overline{b}}^n$ (see (12), (13)), matrix $\tilde{b}$ (see (14) with comments), and scalars $r_m$ , $s_m$ are also dependent variables and subjected to variation. Therefore, the total expression for the nodal force (16) variation is too cumbersome and consists of two items: symmetrical and non-symmetrical. In practice of metal forming analysis, a non-symmetrical item can be rejected without the loss of convergence. A symmetrical positively defined item of the stiffness matrix of the nodal finite element is obtained by the first item in the right hand of (19) and is of the form (p and q are node numbers of complex $C_n$ )

${}^{(n)}K{}_{qk}^{pi}=\frac{{\pi }^2}{3V_n}{\overline{b}}_{pj}^{\overline{n}}(C_{\overline{n}}^{ijkl}-{\stackrel{⌢}{C}}_{\overline{n}}^{ijkl}){\overline{b}}_{ql}^{\overline{n}}(1\le i,j,k,l\le 2)$ (20)

while if $i=1$ or $k=1$ then some special items should be added to the stiffness matrix respectively:

$\Delta {}^{(n)}K{}_{qk}^{p1}=\frac{2{\pi }^2}{9V_n}{\overline{b}}_{ql}^{\overline{n}}(C_{\overline{n}}^{33kl}-{\stackrel{⌢}{C}}_{\overline{n}}^{33kl}){\tilde{b}}_p^{\overline{n}}$ or $\Delta {}^{(n)}K{}_{q1}^{pi}=\frac{2{\pi }^2}{9V_n}{\overline{b}}_{pj}^{\overline{n}}(C_{\overline{n}}^{ij33}-{\stackrel{⌢}{C}}_{\overline{n}}^{ij33}){\tilde{b}}_q^{\overline{n}}$ $(1\le i,j,k,l\le 2)$ (21)

and if $i=k=1$ then one more item should be added:

$\Delta {}^{(n)}K{}_{q1}^{p1}=\frac{4{\pi }^2}{27V_n}{\tilde{b}}_p^{\overline{n}}(C_{\overline{n}}^{3333}-{\stackrel{⌢}{C}}_{\overline{n}}^{3333}){\tilde{b}}_q^{\overline{n}}$ (22)

When assembling the global stiffness matrix, the stabilizing matrices are also taking into account (they are calculated according to classical expressions in Zienkiewicz (1971) and constitutive equation (18) ). This completes the nodal averaged axisymmetric finite element formulation.

3 NUMERICAL EXAMPLES

The following two example problems demonstrate the abilities of the new finite element. The first test is a benchmark involving large strain, elastoplasticity, and specific “stone wall” boundary conditions. The second example shows a possible industrial application when a three dimensional model is reduced to an axisymmetric one with considerable cost saving. As for the stabilizing tensor $\stackrel{⌢}{C}$ , the natural choice $\stackrel{⌢}{C}~C$ is not optimal. Better results may be obtained by letting (see Puso and Solberg (2006) )

${\stackrel{⌢}{C}}^{ijkl}=p(\stackrel{⌢}{\lambda }{\delta }^{ij}{\delta }^{kl}+\stackrel{⌢}{\mu }({\delta }^{ik}{\delta }^{jl}+{\delta }^{il}{\delta }^{jk}))$ (23)

where p is a small penalty parameter (p=0.05 is recommended) and Lame’s constants $\stackrel{⌢}{\lambda }$ and $\stackrel{⌢}{\mu }$ are defined in the following way. For elastic material with Lame’s constants $\lambda$ and $\mu$:

$\stackrel{⌢}{\mu }=\mu$ and $\stackrel{⌢}{\lambda }=\min (\lambda \mathrm{,25}\stackrel{⌢}{\mu })$ (24)

For plastic materials, a shear modulus based on the plastic tangent modulus E _T is used:

$\stackrel{⌢}{\mu }=E_T/2$ and $\stackrel{⌢}{\lambda }=\min (\lambda \mathrm{,25}\stackrel{⌢}{\mu })$ (25)

A contact problem statement conforms to Morrev (2007 , 2009 , and 2011 ). Only the symmetrical positively defined items (20) , (21) , (22) of the finite element matrix are evaluated.

3.1 Cylindrical billet upsetting

In this example, which is adopted in Puso and Solberg (2006) , performance on a large deformation, confined plasticity problem is evaluated. The undeformed steel (E =2.987×10⁴ ksi, ν=0.3) billet is 2 in in height and diameter and quarter symmetry triangular mesh model is used for simulation as seen in Figure 3 . The billet is squashed by two parallel rigid plates to thickness 0.5 in. A contact surface of type “stone wall” (i.e. with infinite frictional coefficient) is used to vertically constrain the outer boundary as it expands (bulges) horizontally. The material is elastoplastic with the power law hardening of the form

Figure 3
Upsetting of cylindrical billet. Effective plastic strain is shown for three different thicknesses: (a) 2 in; (b) 0.82 in; and (c) 0.5 in.

${\sigma }_y({\epsilon }^p)=112{({\epsilon }^p+{\epsilon }^0)}^{0.227}ksi,{\epsilon }^0={(41/112)}^{1/0.227}\approx 0.0113$ (26)

In this example, the tangent modulus E_T , which is used in stabilizing tensor $\stackrel{⌢}{C}$ , was the tangent of the power law equation evaluated at ε_p =0:

$E_T=\frac{d{\sigma }_y}{d{\epsilon }^p}(0)=112(0.227){(0.0113)}^{-0.773}\approx 813ksi$ (27)

The only difference from the original work by Puso and Solberg (2006) in the problem statement is the optimal penalty parameter value p=0.01 instead of p=0.05. That is, the axisymmetric finite element and the tetrahedral one are not identical in spite of identical meshes and common theoretical framework! The results are presented in Figure 3 where the smoothing procedure (see Section 1) was applied twice to both the isoline patterns. It is difficult to compare these isolines with the results by Puso and Solberg (2006) because in that work color pictures were used to visualize the tensor fields. Nevertheless, the qualitative similarity is evident.

3.2 Deep rolling of a steel rod

A surface hardening process of a steel rod by means of deep rolling (DR) is studied. DR is an extremely complex process (see for example Radchenko et al. (2015) and Gryadunov et al. (2015) ) because of high gradient values of tensor fields and due to the induced non-uniformity of a material. As a rule, such problems are solved in a three-dimensional statement (despite a workpiece looks like an axisymmetric body) when a small representative part of a detail is considered together with some reasonable boundary conditions. In contrast to this approach, the proposed axisymmetric finite element permits to consider the entire detail in some events. The point is that a shape of a contact zone depends on ratio R/r where R is longitudinal roller’s radius and r is transversal one. If $R\approx r$ then the contact zone is, roughly speaking, round and the material flows in all directions from under the roller. But, if $R>>r$ then the contact surface is of the form of a very oblong oval hence the material flow is negligible in the direction of rolling and predominates in the direction of the workpiece axis, i.e. the flow character is the same as in axial symmetric models. Consequently, a three dimensional model can be reduced to axial symmetric one if radius R is sufficiently large. The computation time diminishes drastically in such case. Turn to an appropriate example. A steel rod of radius 20 mm and of length $l=100$ mm is rolled by two rollers of transverse radius $r=0.908$ mm and of longitudinal one $R>>r$ while the rollers move in opposite directions ( Figure 4 ). Note once again that if $R\approx r$ then the problem cannot be considered as an axial symmetric one. Rolling depth $h=0.15$ mm; rolling step d=0.5 mm; number of passes is equal to 1. Material is steel AISI 4340; hardening law is $\sigma ={\sigma }_0+b{\epsilon }^n$ (i.e., slightly different from (26), (27)) where ${\sigma }_0=792$ MPa; $b=510$ MPa; $n=0.26$ ; Young’s modulus $E=210$ GPa; Poisson’s ratio $\nu =0.324$ .

Figure 4
Symmetrical deep rolling

The operation shown in Figure 4 is not realistic but this process possesses the bilateral symmetry therefore the total number of nodes may be reduced twice. In practice, the working is carried out by a single roller or by two-three rollers rotating in a common plane. The finite element mesh is presented in Figure 5 where the finest discretisation is marked out by black; the size of this region is 5×1mm and the splitting numbers here are 50×10; the total number of nodes is 2349. The feature of this example is too high gradient values of tensor field (due to local character of deformation) such that a great value $\stackrel{⌢}{\mu }=E/2$ in (23) , (24) , (25) appears to be necessary instead of $\stackrel{⌢}{\mu }=E_T/2$ . The penalty parameter value p=0.05 instead of p=0.01 in the previous example. The residual stress component along the rod axis is to be found.

Figure 5
Finite element mesh

A like problem was investigated in Majzoobi et al. (2010) with the only difference: longitudinal radius was equal to $R=2.5$ mm. The results of both the works are compared in Figure 6 . The peak values are almost equal. The difference between the two plots is explained by different values of radius R. It is obvious that if $R>>r$ then the contact zone is of the form of a very oblong oval and is much larger than in case $R=2.5$ mm. Correspondingly, the penetration depth of plastic deformation is larger too when $R>>r$ . It should be pointed out that the three-dimensional model in Majzoobi et al. (2010) was calculated up to depth 3 mm only and in the present work an extrapolation of those results up to depth 4 mm is used. The developed finite element permits to evaluate the residual stress without any restriction in depth and a proper result is presented in Figure 7 .

Figure 6
Axial residual stress up to deph 4mm

Figure 7
Axial residual stress along the entire deph

4 CONCLUSIONS

A stabilized nodal averaged axisymmetric finite element was developed. It based on the virtual work principle for large irreversible strain together with a nodal averaged approximation and an approximate formula for a ring of triangular profile. An additional stabilizing stress tensor appears to be necessary for incremental numerical analysis. Simple expressions for nodal force and stiffness matrix permit to obtain fast algorithms and short program codes. The finite element can be exploited successfully for numerical solving the most complicated task of solid mechanics – multiple periodic local deformations as it takes place in the DR process. The axisymmetric finite element and the tetrahedral one are not completely identical in spite of identical meshes and common theoretical framework. In some cases (when the restriction in radii $R>>r$ holds), it becomes possible to calculate an entire detail under working rather than its small three-dimensional part with considerable reduction in computation time and input data preparing.

ACKNOWLEGEMENTS

Special thanks to E. Yu. Beregovaya for her assistance in program code debugging and for developing of pre- and post-processors. The work was carried out within the bounds of the Russian Federation government target 1.5265.2017/БЧNº1.5265.2017/8.9.

• http://dx.doi.org/10.1590/1679-78254349

References

• Areias, P. M., Antonio, C. A., Fernandes, A. A. (2003). Analysis of 3d problems using a new enhanced strain hexahedral element. International Journal for Numerical Methods in Engineering. 58: 1637–1682.
• Beissel, S., Belytschko, T. (1996). Nodal integration of the element-free Galerkin method. Computer Methods in Applied Mechanics and Engineering. 139: 49–74.
• Bonet, J., Marriot, M., Hassan, O. (2001). An averaged nodal deformation gradient linear tetrahedral element for large strain explicit dynamic applications. Communications in Numerical Methods in Engineering 17: 551–561
• Brenner, S. C., Scott, L. R. (2008). The Mathematical Theory of Finite Element Methods (Third Edition). Springer. New York.
• Broccardo M., Micheloni M., Krysl P. (2009) Assumed-deformation gradient finite elements with nodal integration for nearly incompressible large deformation analysis. International Journal for Numerical Methods in Engineering. 78:1113-1134.
• Chen, J. S., Wu, C. T., Yoon, S., You, Y. (2001). A stabilized conforming nodal integration for Galerkin mesh-free methods. International Journal for Numerical Methods in Engineering. 50: 435–466.
• Cueto, E., Chinesta, F. (2015). Meshless methods for the simulation of material forming. International Journal of Material Forming. 8: 25–43.
• Dohrmann, C. R., Heinstein, M. W., Jung, J., Key, S. W., Witkowski, W. R. (2000). Node-based uniform strain elements for three-node triangular and four-node tetrahedral meshes. International Journal for Numerical Methods in Engineering 47: 1549–1568.
• Gee, M.W., Dohrmann, C.R., Key, S.W., Wall, W.A., (2009). A uniform nodal strain tetrahedron with isochoric stabilization, International Journal for Numerical Methods in Engineering. 78: 429–443.
• Gil, A. J., Lee, C. H., Bonet, J., Aguirre, M. (2014). A stabilised Petrov–Galerkin formulation for linear tetrahedral elements in compressible, nearly incompressible and truly incompressible fast dynamics. Computer Methods in Applied Mechanics and Engineering. 276: 659–690.
• Gryadunov, I. M., Radchenko, S. Yu., Dorokhov, D. O., Morrev, P. G. (2015) Deep Hardening of Inner CylindricalSurface by Periodic Deep Rolling - Burnishing Process. Modern Applied Science 9: 251–258.
• Lee, C. H., Gil, A. J., Bonet, J. (2014). Development of a stabilised Petrov–Galerkin formulation for conservation laws in Lagrangian fast solid dynamics. Computer Methods in Applied Mechanics and Engineering. 268: 40–64.
• Majzoobi, G.H., Teimoorial Motlagh S., Amiri, A. (2010). Numerical simulation of residual stress induced by roll-peening Transactions of The Indian Institute of Metals 63: 499–504.
• Morrev, P. G. (2007). A Version of Finite Element Method for Frictional Contact Problems. Mechanics of Solids 4: 640–651.
• Morrev, P. G. (2009). A rate variational principle of quasistatic equilibrium for absolutely rigid body in contact problems. Fundamental and applied problems of technique and technology 6: 30–32. [in Russian]
• Morrev, P. G. (2011) A variational statement of quasistatic “rigid-deformable” contact problems at large strain involving generalized forces and friction. Acta Mechanica 222: 115–130.
• Pires, F., de Souza Neto, E. A., de la Cuesta Padilla, J. L. (2004). An assessment of the average nodal volume formulation for the analysis of nearly incompressible solids under finite strains, Communications in Numerical Methods in Engineering. 20: 569–583.
• Puso, M. A. (2000). A highly efficient enhanced assumed strain physically stabilized hexahedral element. International Journal for Numerical Methods in Engineering. 49: 1029–1064.
• Puso, M. A., Solberg, J. (2006). A formulation and analysis of a stabilized nodally integrated tetrahedral. International Journal for Numerical Methods in Engineering 67: 841–867.
• Puso, M. A., Chen, J. S., Zywicz, E., Elmer, W. (2008). Meshfree and finite element nodal integration methods. International Journal for Numerical Methods in Engineering 74: 416–446.
• Quak, W., van den Boogaard, A. H., Gonsales, D., Cueto, E. (2011). A comparative study on the performance of meshless approximations and their integration. Computational Mechanics. 48: 121–137.
• Radchenko, S.Yu., Dorokhov, D.O., Gryadunov, I. M. (2015). The volumetric surface hardening of hollow axisymmetric parts by roll stamping method. Journal of Chemical Technology and Metallurgy 50: 104–112.
• Reese, S. (2002). On the equivalence of mixed element formulations and the concept of reduced integration in large deformation problems. International Journal of Nonlinear Sciences and Numerical Simulation 3: 1–33.
• Reese, S. (2005). On a physically stabilized one point finite element formulation for three-dimensional finite elastoplasticity. Computer Methods in Applied Mechanics and Engineering 194: 4685–4715.
• Simo, J. C, Armero, F. (1992). Geometrically nonlinear enhanced strain mixed methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering. 33: 1413–1449.
• Simo, J. C, Armero, F., Taylor, R. L. (1993). Improved versions of assumed enhanced strain tri-linear elements for 3d-finite deformation problems. Computer Methods in Applied Mechanics and Engineering. 110: 359–386.
• Simo, J. C, Taylor, R. L., Pister, K. S. (1985). Variational and projection methods for the volume constraint in finite deformation elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 51: 177–208.
• Sussman, T., Bathe, K. J. (1987). A finite-element formulation for nonlinear incompressible elastic and inelastic analysis. Computers and Structures 26: 357–409.
• Zhang, Z., Dolbow, J. E. (2017). Remeshing Strategies for Large Deformation Problems with Frictional Contact and Nearly Incompressible Materials. International Journal for Numerical Methods in Engineering. 109: 1289–1314.
• Zienkiewicz, O. C. (1971). The Finite Element Method in Engineering Science, McGraw-Hill, London.

Publication Dates

• Publication in this collection
  2018

History

• Received
  01 Aug 2017
• Reviewed
  18 Feb 2018
• Accepted
  22 Feb 2018