Stochastic Elasto-Plastic Analysis of Necking Bar with Random Tvergaard Coefficients

This paper reports on the computational modelling of static extension tests of the round steel bar. The main objective was to apply the generalised stochastic perturbation technique implemented as the Stochastic Finite Element Method to carry out the numerical simulation of its elasto-plastic behaviour. This approach was based on: the general order Taylor expansion of all input random variables and the resulting state functions of their average means, as well as on the Least Squares Method employed to determine analytical functions of in-between design parameters and the given structural responses. Tvergaard coefficients were assumed as the uncorrelated Gaussian random variables to check the effect of material porosity uncertainty on the statistical scattering of its deformations and stresses. The computational implementation employed the FEM system ABAQUS and computer algebra system MAPLE, including polynomial and non-polynomial local response functions of the displacements, plastic strains and reduced stresses. Moreover, 4-node axisymmetric, continuum, reduced-integration FEM elements (CAX4R) were used in the conducted analyses. The basic probabilistic characteristics of the structural response (expectations, coefficients of variation, skewness and kurtosis) were determined throughout the entire deformation process as the functions of input uncertainty level. The obtained results were finally contrasted with the classical Monte-Carlo Simulation scheme and the semi-analytical technique for input coefficient of variation of porous plasticity coefficients not larger than 0.20.


Introduction
Currently used Eurocode standards allow for using the plastic reserve of the capacity for steel structures. If we take into account the influence of the microstructural damage on the strength of material we cannot use well known Huber-Mises-Hencky strength hypothesis. Material continuum principal does not work here. This occurs when we carry out extension analysis where necking appears. In such a case, we have to use other, more developed damage models allowing defining the material damage range. These material properties can be certainly assumed as a random feature. Determination of such randomness is always challenging, all the more that there exists no scientific output considering such an analysis with random Tvergaard coefficients. Three well-known simulation techniques used in the study were: the semi-analytical method (AM), Monte Carlo Simulation (MSC) and 10 th order perturbation technique (PM) which base on Taylor expansion series [1]. The first four probabilistic moments and coefficients were determined to obtain information about probability density function. The specimens subjected to elasto-plastic test were made of porous cylindrical steel. Tvergaard coefficients i q , as the uncorrelated Gaussian random variables, were to verify the effect of material porosity uncertainty on deformations and stresses. This randomness was performed independently for 1 q and 2 q . Our future scientific plans involve calibration of the presented stochastic numerical model of structural steel containing microdefects by means of full-scale tests. This approach is based on the well-known generalised perturbation method applied together with constitutive Gurson-Tvergaard-Needleman (GTN) material model implemented in ABAQUS Standard system.

Theoretical background
Random variable b can be described by its probability density function   ; the expected value of this variable is defined as [2] where 0 b stands for the average value of the variable b. mth central probabilistic moment of b is given as follows The main idea of the stochastic linearised perturbation approach is to expand all input random variables and all the resulting state functions (displacements, stresses, etc) of the given boundary initial problem via Taylor series about their spatial expectations using the perturbation parameter 1   .
The expected values of horizontal displacements defined by its expansion via Taylor series are: Furthermore, the classical integral Nth central probabilistic moment of the displacement state function x u is given by: 10 th order approach was used, and therefore orders higher than 10 th are neglected. Based on the classical definition of the variance we can determine the coefficient of variation from: Skewness and kurtosis were calculated to determine the type of asymmetry and concentration probabilistic distribution of the state function. Skewness is the ratio of the third central probabilistic moment and the third power of standard deviation as: (6) and kurtosis is given by: It should be mentioned that at this stage the proposed procedure is independent of a choice of the initial probability distribution function. These random characteristics were calculated by stochastic perturbation, Monte Carlo and the semi-analytical method in MAPLE. 100.000 samples were used for MCS with the use of a random number and statistical estimators generator implemented in MAPLE. Since the semi-analytical method uses the same polynomial response functions as the perturbation technique, therefore the probabilistic moments are calculated by analytical integration from definition also implemented in MAPLE. The entire multi-method approach is based on the solution of the iterative deterministic FEM [3] equation: (8) where K denotes stiffness matrix, R is the matrix of nodal loads, q is the solution displacement vector, while  is the tests number necessary for the Response Function Method recovery. The solution of Eq. (8) is provided in ABAQUS system for eleven discrete realisations of each random variable uniformly partitioning the interval ] ; [ . Subsequently, the determination of the probabilistic moments proceeds according to Eqs. (3÷7) owing to the symbolic derivation of all partial derivatives with respect to the given random variable .
b Using computer algebra system MAPLE in all these calculations enabled parameterisation of the resulting moments and probability coefficients as a function of the input coefficient of variation of the random variable. It is helpful in the visualisation of obtained results and implementation of the 10 th order perturbation method. The yield condition of the modified Gurson model, which takes into account Tvergaard coefficients, is given by: where e  is the effective stress according to Huber-Mises-Hencky hypothesis, 0  is yield point, m  is mean stress, while f indicates material voids volume fraction. The constants i q were originally introduced by Tvergaard [4,5] to obtain better agreement with the localised shear band bifurcations than the ones predicted by initial Gurson model. The best predictions of the voids growth and coalescence were obtained . Therefore, the main goal of this study is to determine the effect of these uncertainty coefficients on the behaviour of the tensioned steel specimen.

Numerical simulation
Numerical simulation of static extension tests was carried out for only a quarter of the steel bar due to its axial symmetry (Fig. 1). Geometrical dimensions of the bar were mm R 10 0  and mm L 40 0  which means that ratio was 4 / 0 0  R L were taken according to [6]. A small notch at the bottom surface was performed to ensure necking in the middle of the initial bar; the dimensions of this cut are . 0 Fig. 1 shows a free body cut through the specimen to show its symmetry. The simulations were carried out for the typical S235JR steel, which is one of the basic materials employed in structural engineering, belonging to the carbon steel group, with the maximum content of carbon below %. 2 . 0 The specimen is divided into two parts. The bottom parts, where we expect appearing of the necking, are 400 CAX3 FEM elements, -3-noded linear axisymmetric triangles. The edge of the triangle is mm 5 . 0 ; 1400 CAX4R FEM elements were used for the rest of the specimen [7], which are standard 4-noded elements with reduced integration (Fig. 2). The two types of discretisation were implemented on account of the fact that we expected considerable displacement and stress distribution in the region where the notch is located. The implemented methods should improve our calculations and provide more extensive information on this phenomenon. Due to the symmetry of the specimen natural kinematic boundary conditions are 0     10 polynomial response functions were recovered to describe the effect of the Tvergaard coefficient i q b  on this structure behaviour. Horizontal maximum displacement at point A was chosen as the state function (Fig. 1) (Fig. 4).  25] with expectation equal to 0 .
; the displacements were written after each 10% of the analysis progress. These polynomial response functions (order from the 2 nd to the 10 th ), which describe the sensitivity of the specimen with respect to the assumed random variables, were obtained by means of Taylor-Newton-Gauss algorithm [9] (TNGA). This method is helpful in solving non-linear problems and was implemented in the computer algebra system MAPLE. Figs 6÷13 show expectations, coefficients of variation, skewness and kurtosis of the extreme x u displacements at point A (Fig. 1) of the specimen.   (Fig. 9) are close to zero as well and rapidly increase at the end of the analysis, where the deformations are substantial. Three probabilistic methods do not coincide at the end of the analysis when the neck is the biggest and numerical analysis is highly nonlinear. Skewness (Fig. 10) is close to zero for all three methods if we narrow the input coefficient of variation of random variable     Considering skewness distribution (Fig. 11), one can see fluctuations at the end of the analysis. We notice convergence of all three methods up to 80% of analysis progress. It follows that the displacements can be recognised as a symmetric distribution. After that moment skewness takes positive and negative values and cannot recognise asymmetry of x u displacement distribution. Most values of the displacements are concentrated around their means because their kurtosis (Fig. 12) takes positive values throughout the entire analysis.
It needs to be emphasised that enormous fluctuations of kurtosis, which takes extreme two times for analysis progress, are equal to 10% and at the end of computations were obtained using the semi-analytical method (AM).  The semi-analytical method (AM) and Monte-Carlo method (MCS) produce extreme values at the end of the numerical analysis. Skewness and kurtosis (Fig. 13) increase enormously at the end of the entire implementation process, unlike their increase as the function of input uncertainty level in the elasto-static analyses.

Conclusion
The paper demonstrated a four-moment stochastic numerical analysis with elasto-plastic material exhibiting structural defects, conducted with the generalised stochastic perturbation approach. The computations were carried out by means of common FEM software ABAQUS and computer algebra system MAPLE, which are crucial to Stochastic Finite Element Method (SFEM). Different nature response functions had to be employed to efficiently fit the response of the structure (Figs 4, 5) to model the necking problem with porous metal plasticity. As it was demonstrated, the computer simulation based on the Taylor series perturbation approach (PM) can be an efficient alternative for the well-known Monte-Carlo Simulation (MCS). Convergent results based on three numerical methods were obtained at the beginning of the analysis when elastic deformations appear. At the end of the analysis, when the neck occurs, MCS and the analytical method produced different results than the perturbation method. It was observed that a bigger neck appears when 1 q is assumed as the input Gaussian random variable and this conclusion refers to all three stochastic computer methods. The analytical description of the displacement distribution x u obtained in this study can be further used in the stochastic reliability analysis of the carbon steel structures with porosity implemented according to the Gurson-Tvergaard-Needleman material model for both First and Second Reliability Methods. However, the very interesting case would be random porosity of the defects f, which would be further justified by experimental measurements and technological processes. The key disadvantage of this model consists in that the same porosity level may be obtained with various combinations of defect size and number, which usually leads to various physical phenomena.