Configurational-force-based finite element mesh refinement for elastic-plastic problems

In this paper the computation of configurational forces in case of elastic-ideally plastic material will be examined. Numerical computation of the error in configurational forces will also be introduced in elastic and plastic domain. It will be shown that the so-called r-adaptive mesh refinement procedure [1] is also applicable for elastic-plastic problems as well as the configurational force driven h-adaptive scheme. In some special examples configurational forces are computable in analytic way. This is useful to compare the solution with numerical results, therefore validating the finite element procedure. Two plane problems will be considered where analytical solutions are known. The first one is the thick walled tube model loaded by internal pressure. Second one is an artificial problem where the displacement field assumed to be known in every point of the domain considered. According to the papers Krieg [5] and Szabó [8] analytical solution is obtainable for stress and strain distributions, if the time derivative of the strain is constant. R− and h−adaptive procedures will demonstrated on these two examples.


Introduction
Configurational mechanics presented by Eshelby in the early fifties is widely used on several areas of theoretical continuum mechanics and computational mechanics as well.Usage of material-or configurational forces in FEA is a very popular topic.One of them is the so called r-adaptive FE mesh refinement strategy which was introduced by Braun (1997) [1].According to definition the configurational force [4] should be zero in the interior of the material when the material is assumed to be homogeneous and isotropic [6].However, in case of finite element computation this requirement won't be fulfilled and nodal configurational forces appear on interior nodes too.The origin of these forces is the fact that the results of finite element solution depends on the initial nodal configuration.In case of optimal mesh internal configurational forces are vanishing and the total potential of the structure will be minimal.Considering that configurational force vectors point to direction of increasing total potential, moving the nodes in opposite direction the optimal mesh is obtainable.Since configurational force indicates the error on finite element mesh, computation of this quantity is also suitable to drive h-adaptive mesh refinement.These methods were demonstrated in several articles for elastic case, e.g.[6], [1], [3].In this paper will be shown that they are applicable also in case of elastic-ideally plastic deformation.Small strain and linear elastic behaviour will be considered.

Theoretical background
Brief conception of configurational force -sometimes called material force -could describe as an energy change in the system considered according to the structure of it.

Elastic domain
Strain energy density for linear elastic, isotropic and homogeneous material could have the form (W(ε) where ε is the elastic strain.The gradient of this function is which after rearrangement [6] follows to the so-called material equilibrium equation where is the Eshelby tensor.Configurational force is obtainable with integration of (Eq.2) over domain Ω and this should be zero under circumstances mentioned above [6].Therefore, if the strain energy density has the form then configurational force is

Plastic domain
The strain energy density of an elastic-plastic material could be written as the sum of an elastic and a plastic part according to [7] as The elastic part is function of ε e elastic strain (homogeneous and isotropic material) and the plastic part is function of ε p equivalent plastic strain, respectively.In the following parts ideally plastic behaviour is considered therefore the second part will vanish.Total strain tensor could be decomposed into an elastic and plastic part as and for ideally plastic materials the second term doesn't exist.
Substituted (6) into (7) and taken into account that σ = ∂Φ e ∂ε e after rearrangement div (Wδ − σε) +σ : ∇ε p g = 0. ( Integrating this equilibrium equation over the plastic domain Ω the sum of the two parts should be zero.Because of the numerical error of finite element approximation the equation won't be fulfilled therefore the configurational force produced by this error is where and Considering that for elastic case the total strain ε will be equal to the elastic strain ε e the above equation leads to (4) in the elastic domain.

Configurational force driven FE mesh refinement
First step to compute nodal configurational forces is to evaluate first and second part of expression (8).The sum of them is the error configurational volume force g error .In view of nodal values of this quantity discrete nodal configurational forces obtainable through integration with use of Gaussian quadrature.
After the computation is done for every elements the global configurational force is necessary to assembly from element values.Therefore the global error configurational force vector will be Using the nodal configurational forces the r-adaptive refinement will be driven according to expression where c is an appropriately selected parameter which drives the changes in nodal positions and X contains the nodal coordinates.The value of c could be a small fraction of the quotient of maximal configurational force norm and minimal element size.This expression is only applicable on free degrees of freedom which means that boundaries of the model have to remain intact during the mesh reconfiguration.In case of hadaptive method the local mesh refinement is also driven by the norm of nodal configurational force.If the value of this norm exceed a required limit than it is necessary to refine the mesh around the nodes in question.

Numerical examples 4.1 Thick walled tube -r-adaptivity
Analytic solution for elastic-plastic deformation of thick walled tube loaded by internal pressure is well-known.The model and the material parameters are shown in Fig. 1.Elasticideally plastic material, plane strain and Tresca yield criterion considered.Stress and strain distribution [2] in the plastic zone in polar coordinate system are where the displacement function relating to the plastic domain is In (17) G denotes the shear modulus.According to these solutions the volume configurational force by analytical evaluation of the two parts of (8) will be This is the radial component.Others will be zero because of the symmetry conditions.The diagram of this function is shown in  strategy in a real plane problem it is necessary to find a method to generate analytic solution for stress and strain field which depends on two space variables.If the displacement field is known in every point of a given domain and this displacement is a linear function of time then it is possible to get the analytical solution according to [8].However, such "artificial" mechanical models probably have nothing to do with real applications but they could be appropriate to compare with solutions of numerical methods.Furthermore, it is possible to create solutions using the Mises yield criterion and also hardening effect could be taken into account.The following expressions are written according to [8].Displacement functions in x and y direction let u(x, y, t) and v(x, y, t).If the displacement vector is u = (u(x, y, t), v(x, y, t), 0) the total strain tensor in linear case has the form Hooke's law establish the relation between stress and strain with the formula.The deviatoric part in its eigensystem is and the Mises equivalent stress as the function of deviatoric stress tensor will be Strain deviator has the same form as equation ( 21), so and this tensor also considered in its eigen coordinate system.The yield stress from an arbitrary, purely elastic state is computable as In the above equation k∆e is the elastic part of the strain deviator.
For the computation of k we may use the formula (by virtue of Prandtl-Reuss plasticity theory).As further deformation occurs the stress deviator could be written as the linear combination of yield stress and the deviatoric strain increment, where the two coefficients are In equations (27) R 0 is the yield stress and for ideally plastic material R 0 = R n+1 .Furthermore, angles ψ 0 and ψ n+1 are and respectively.Substituting the angles into equations ( 27) and the parameters from ( 27) into (26) we get the stress deviator inside the plastic domain.The stress tensor in our original coordinate system will be Fig 4 shows a specific example where the shape of the plastic domain would be determined by the method described above.Material model assumed to be linear elastic and ideally plastic with Mises yield criterion.Plane strain is considered.Displacement field would be prescribed as u(x, y, t) = t sin x • cos y v(x, y, t) = t sin y • cos x .
Material parameters and the model shown in Fig. 4a with the plastic zone at t = 0.0008.Fig. 4b shows the nodal error configurational force distribution at a given path.Increasing the mesh density the error will be decreasing.In this case the mesh refinement was uniform but the application of local refinement could increase the effectiveness of the computation.In this case only certain elements are affected on which the error exceed a prescribed limit.

Conclusions
R− and h−adaptive finite element mesh refinement methods were proposed which are based on numerical configurational force computation.Brief conception of these two methods were introduced associated with two different elastic-plastic plane problem which also has analytical solution for stress and strain distribution, respectively.The so called r−adaptive scheme is also applicable for elastic-plastic analysis.The computation of nodal error configurational force is slightly different as in purely elastic case.However, application of this method require more computation time because the optimal nodal configuration only achievable through iteration and it is necessary to solve the original problem in each step.In case of elastic plastic analysis this  happens through several increment.More promising approach to use configurational forces to indicate error for h−adaptive scheme.In this case the dimension of the problem will be greater than by r−adaptivity but only a few iteration step could be enough to achieve the required accuracy.

Fig. 2
compared to the numerical solution.The difference between the two quantity indicates the error in nodal configuration of the mesh.Finite element solution was carried out with use of quadratic and axisymmetric quadrilateral elements.Discrete error configurational forces are shown in Fig.3.After 35 iteration step with 10 elements along radial direction of the wall the error significantly decreased in the plastic zone as it would be used double element density.

Fig. 4 .
Fig. 4. a.) Plate model with material parameters and prescribed boundary conditions u(x, y, t) = t sin x • cos y and v(x, y, t) = t sin y • cos x b.) Nodal configurational error force along y = π 4 line