Experimental and Computational Study of Ductile Fracture in Small Punch Tests

A unified experimental-computational study on ductile fracture initiation and propagation during small punch testing is presented. Tests are carried out at room temperature with unnotched disks of different thicknesses where large-scale yielding prevails. In thinner specimens, the fracture occurs with severe necking under membrane tension, whereas for thicker ones a through thickness shearing mode prevails changing the crack orientation relative to the loading direction. Computational studies involve finite element simulations using a shear modified Gurson-Tvergaard-Needleman porous plasticity model with an integral-type nonlocal formulation. The predicted punch load-displacement curves and deformed profiles are in good agreement with the experimental results.


Introduction
Small punch (SP) testing is used at high homologous temperatures for investigation of creep properties such as rupture time and minimum creep rate or at low homologous temperatures for investigation of other mechanical properties such as yield stress, ultimate stress or fracture toughness. The test requires smaller specimen sizes as compared to standard mechanical tests. Hence, it allows investigations of regions with gradients properly, such as heat affected zones [1,2]. Due to its small size the test is almost non-destructive. Thus, it eliminates the need for repairing the component after sample removal, which is another advantage of the SP test compared to the standard tests. During an SP test the sample deforms by different deformation mechanisms initially by elastic and plastic bending and followed by membrane stretching. Hence, a rather complex multiaxial stress state, evolving with puncher displacement, occurs inside the sample. The sample thickness closely affects the occurrence of these deformation modes. If the thickness is increased then shearing mode becomes more dominant as compared to membrane stretching [3]. In a recent work [4], the authors applied Gurson-Tvergaard-Needleman/Ritchie-Knott-Rice (GTN/RKR) approach to plain and notched SP disks for a wide temperature range. The current study aims at a detailed investigation of deformation mechanisms as well as crack initiation and propagation at room temperature SP tests of P91 steel with different disk thicknesses. To this end, the problem is studied by experimental and numerical analyses in order to: • Further exploit multiaxial SP testing for numerical model validation, • Consolidate the parameter identification through application of the model to a wide range of disk thicknesses, • Elucidate the effects of shear and initial porosity on damage progression and failure, • Derive robust predictions for the yield stress (YS) and ultimate tensile stress (UTS) for future applications.
On the experimental side, SP tests are carried out using disks of P91 material which is a heat resistant steel widely used in power plants. Due to this reason it is frequently preferred as the test material in the SP method [1,2,5,6]. The SP disks are produced with varying thicknesses from 0.2 mm to 1 mm were tested up to complete fracture. At all thicknesses, due to lack of notch, large scale yielding conditions prevail (at least before intense localization with voidage) and the observed consequent fracture mode is ductile. For smaller thicknesses the fracture is led by diffuse necking followed by an intense localization throughout the section. For larger thicknesses, on the other hand, an initial diffuse necking pattern is not observed since shear deformation prevails rather than membrane stretching. Unlike a simple tension test where the cracks start inside the specimen, the crack initiates from the bottom surface at a distance where thinning takes place. As will be shown, with increasing thickness the fracture occurs at larger punch displacements regarding both the computational and experimental curves. Scanning electron microscope (SEM) studies of the fracture surfaces show signs of a mixed mode I (normal) mode II (shear) fracture with relatively flat dimple walls elongated along the fracture surface. Fracture is primarily due to voidage nucleated at M 23 C 6 carbides and/or MX precipitates coalesced under severe shearing.
On the computational side, finite element simulations are carried out for the given thickness range. To this end a shear modified GTN model with strain hardening and strain rate hardening is implemented to model the ductile fracture based on void nucleation, growth and coalescence. The origin of the model goes back to Gurson [7] and later modifications were introduced by Tvergaard and Needleman to better account for cavity growth [8][9][10]. More recently, a shear modification has been proposed in [11] to accommodate softening by inter-void linking under low stress triaxiality conditions which is not accounted for in the original formulation. Pathological mesh sensitivity pertaining to softening is remedied by an integral-type nonlocal formulation applied to each additive void volume fraction rate component similar to [12]. This formulation requires a characteristic length hence allowing incorporating the size effect in the model. This model is implemented as a VUMAT user defined material subroutine for ABAQUS/EXPLICIT. Prior to applications to the SP test, the implementation is validated using single element tests under dilatation and simple shear conditions for which analytical solutions exist. The effectiveness of the implemented delocalization scheme is also demonstrated by plane strain tensile tests with different mesh sizes. These studies are followed by SP simulations with 2D axisymmetric models. Model parameters for P91 steel including the associated length scale are identified using quantitative measurements, i.e., metallographic images as well as inverse methods relying on experimentally determined force-displacement curves. A good agreement of not only the force-displacement curves but also the deformed profiles of the failed samples in between the computationally and experimentally determined results is obtained. A parameter sensitivity analysis is also conducted where the influence of shear damage parameter k w and the initial porosity f 0 are investigated. Since the gradient of solution dependent fields are relatively low during the tests in these investigations, the effect of nonlocal regularisation is not observed until the onset of localization. It is believed that the nonlocal modeling framework developed will be helpful to further analyze the effects of strain rate, specimen size, puncher head geometries and especially notched specimens where higher gradients of field variables could prevail.

A Word on Notation
Assuming a, b and c as three second order tensors, together with the Einstein's summation convention on repeated indices, c = a · b represents the single contraction product with c ik = a ij b jk . d = a : b represents the double contraction product with d = a ij b ij , where d is a scalar. E = a ⊗ b, F = a ⊕ b, and G = a b represent the tensor products with E ijkl = a ij b kl , F ijkl = a ik b jl , and G ijkl = a il b jk , where E, F, and G represent fourth order tensors. a and a −1 denote the transpose and the inverse of a, respectively. ∂ a b denote ∂b/∂a. dev (a) and tr (a) stand for the deviatoric part of and trace of a, respectively, where dev (a) = a − [1/3]tr(a) 1, with tr(a) = a ii and 1 denoting the identity tensor. sym (a) and skw (a) denote symmetric and skew-symmetric portions of a.ȧ gives the material time derivative of a. x = [1/2] [x + |x|]. Finally, a gives a represented at the rotationally neutralized configuration.

Fundamental Kinematical Assumptions and Hypoelastic Plasticity
Let F := ∂ X x define the deformation gradient of the nonlinear map ϕ : Ω 0 × R → R 3 where X and x := ϕ(X, t) denote the particle positions at the reference (undeformed) configuration Ω 0 and current (deformed) configuration Ω respectively. Then, l :=Ḟ · F −1 = ∂ x v denotes the spatial velocity gradient, with v =ẋ. An additive split of d := sym (l) into elastic and plastic parts (Various forms of rate additive splits can be recovered through a multiplicative decomposition of F into elastic F e and plastic F p parts with F := F e · F p making use ofḞ =Ḟ e · F p + F e ·Ḟ p and For small elastic strains and rigid body rotations with F e 1 one governs l p → L p . Using F e := R e · U e and assuming small elastic strains but finite rotations, i.e., U e 1, supplies l p → R e · L p · [R e ] −1 . Taking the symmetric part of both sides one reaches d = d e + d p . A small elastic strain assumption makes a hypoelastic constitutive equation in terms of d e plausible.) is hypothesised to reach d = d e + d p . We then introduce the rotationally neutralized rate of deformation tensor˙ defined aṡ with˙ e := R · d e · R and˙ p := R · d p · R. Here, R is the rotation tensor found through the polar decomposition of the deformation gradient with F = R · U. Similarly, a pull back operation on the Cauchy stress tensor σ with the rotation tensor gives its rotationally neutralized counterpart viz σ := R · σ · R whose material time derivative˙ σ is postulated to obey the following hypoelastic relatioṅ with C e : where K and µ denote the bulk and the shear modulus, respectively.

Shear Modified GTN Porous Plasticity-Local Formulation
The model used here is Gurson's dilatant plasticity model [7], which is extended by parameters q 1 , q 2 and q 3 in [8,9] to account for a better agreement with the computational analyses of various void distributions, and by the bilinear function f * ( f ) in [10] to account for rapid void coalescence prior to failure. The hydrostatic stress dependent flow potential Φ p is then formulated as with where σ eq is the (macroscopic) equivalent von Mises stress, σ y is the (matrix) yield stress, σ m is mean stress, f is the void volume fraction and q 1 , q 2 and q 3 are material parameters. f c denotes the critical void volume fraction at incipient coalescence, f f the fraction at failure and f * u = 1/q 1 . The viscoplastic hardening of the material matrix is described by σ y which accounts for strain and strain rate dependence. Hence, letting e p denote the equivalent plastic strain andė p its rate, we assume the following multiplicative form where h y and r y respectively denote the functions of strain hardening and strain rate hardening where σ y0 , σ y1 , σ y∞ , h 0 , h 1 , m, n and e p 0 are plastic strain hardening parameters. C andė p 0 , on the other hand, are parameters controlling the rate dependence.
Associated plastic flow rule gives the plastic rate of deformation tensor at the rotationally neutralized configuration as˙ whereγ denotes the plastic multiplier. The equivalent plastic strain rate, using the plastic work equivalence with [1 − f ] σ yė p = σ :˙ p , readṡ The void volume fraction evolution involves nucleation and growth. The rate of the total void volume fraction is formulated additively in terms of void nucleation rateḟ n and void growth rateḟ ġ f =ḟ n +ḟ g .
Forḟ n we assume a strain dependent void nucleation [13] where e p N and S N respectively denote the mean equivalent plastic strain at the incipient nucleation and its standard deviation. f N denotes the total source for void volume fraction of nucleation.ḟ g is further split into two parts [11] whereḟ g normal accounts for the void growth under hydrostatic stresses whose formulation simply uses the mass conservation viz.ḟ whereasḟ g shear accounts for softening effects associated with void distortion and void interaction with material rotation under shear stress states wherė Here, k w is a material parameter which scales shear damage rate with a suggested interval with 0 ≤ w ≤ 1 distinguishes the states of axisymmetric stress from those of generalized shear on the Π−plane. Here, w = 0 through eq for all axisymmetric stress states, whereas w = 1 through J 3 = 0 for the states of generalized shear. This modification due to [11] is motivated by the reported experimental evidence for low triaxiality fracture development in, e.g., [14,15], for which the original GTN model falls short in predictive capability.
The fracture responses for the shear modified Gurson model are given in Figure 1 for linear and uniform strain paths under plane stress state considering linear isotropic hardening plasticity with hypothetical parameters. As the figure clearly depicts the monotonic dependence of the fracture strain on the stress triaxiality is suppressed for k w = 0. Further, there occurs considerable shrinkage in the admissible range of deformation for generalized plane strain paths, i.e., pure shear and plane strain loading paths with vanishing stress triaxiality ratio σ m /σ y → 0. On the other hand, strain paths associated with the axisymmetric stress states, i.e., uniaxial and biaxial loading paths respectively with σ m /σ y = 1/3 and σ m /σ y = 2/3, are not affected by the shear correction since for these paths the shear fracture controlling parameter w(dev σ) becomes zero with vanishing third invariant of the deviatoric stress tensor: J 3 = 0. Applications of the model to formability exhaustion of DC04 steels during sheet-bulk forming and to the problem of severe plastic localization bands initiated at free surfaces during free bending are given in [16,17], respectively.

Integral-Type Nonlocal Extension
The motivation behind integral-type nonlocal formulations is twofold. On the mechanical side, with micro-void and micro-crack interactions being their micromechanical motivation, integral-type nonlocal formulations constitute distributed damage models capable of reproducing size effects. In view of finite element analysis with damage models, they remedy the pathological mesh dependence of the local solution where the size of the process zone and associated energy dissipation per unit volume is dictated by the discretization. With this motivation, in the current study, an integral-type nonlocal formulation is adopted which relies on the following delocalization operation Here y represents the location vector and V the volume at the current coordinates. If ω (x, y) denotes the bell-shaped nonlocal weight function the normalized weight function ω (x, y), which remedies any inconsistency pertaining to the unrestricted averaging domains extending over the problem boundary, reads As long as the boundaries are not violated V ω (x, y) dV (y) is a constant. R in Equation (15) denotes the interaction radius which constitutes the characteristic length. For R → 0 a local formulation is recovered. In practical applications, R is related to the material microstructure, e.g., four times void size or half void spacing for ductile fracture mechanism [18]. Delocalization could be applied either directly to the kinematic fields; see, e.g., [19], or to the damage variables which control softening, see, e.g., [12,18,20,21]. In the current study, delocalization of each additive rate component of the void volume fraction is applied. Hence, v in (14) is substituted by the damage rate component as followṡ which finally adds up to the total nonlocal porosity rateḟ nonlocal , viz.
One notes that for spatially uniform porosity rate distributions delocalization is not effective. Hence, local and nonlocal theories differ only if field gradients exist. The material model is implemented as a VUMAT material subroutine for ABAQUS/EXPLICIT. The algorithmic details are given in the supplementary materials. The natural setting of VUMAT alone is not appropriate for nonlocal formulations and it allows only local effects to be modeled. To this end, in our treatment we have additionally implemented a VUSDFLD routine.

Material Employed
The material used in this study is modified P91 which has a tempered martensitic microstructure. An optical image revealing the microstructure is given in Figure 2a

Experimental Setup
SP tests were performed using the experimental setup in Figure 2b. The major components of the set up are a specimen holder to ensure the tight clamping of the SP disk, a puncher with a hemispherical head of 2.5 mm diameter for the central loading of the disk and two Linear Variable Displacement Transducers to measure the displacement of the puncher. The SP specimens were circumferentially clamped to prevent slippage of the SP specimen during the test. The aperture of the receiving die was 4 mm in diameter with an edge chamfer of 0.2 mm. The crosshead displacement rate was 0.005 mm/s. All the tests were carried out at room temperature (25 • C) in accordance with the Code of Practice for Small Punch Testing [22]. The typical output from an SP test is a force-displacement curve.

Crack Propagation and Fracture Surfaces
During SP testing at room temperature, the fracture is ductile and proceeded by uniform necking. The crack initiates from the bottom surface at a distance where necking takes place. Then, it propagates in direction of the maximum equivalent plastic strain through the thickness and follows a circumferential path along this necking region [23,24]. This is valid regardless of disk thickness.
For most metals, voids nucleate from inclusions and secondary phase particles by particle cracking or interface decohesion with increasing plastic strain. If particles are not large like MnS inclusions, the voids nucleate by debonding of the particle-matrix interface and grow with the plastic deformation in the matrix. The resulting fracture surface exhibits a dimpled structure with many microvoids. In Figure 2c the fracture surface of the SP disk with 0.5 mm thickness is presented. This image reveals a dimpled fracture surface consisting of high density of small microvoids and lower density of relatively large and deep ones. The approximate distance between the large dimples is found to be 30-35 µm whereas the distance was 2-5 µm for the small dimples. The initiation of small microvoids is attributed to MX precipitates distributed in the matrix which are higher in density. The larger microvoids presumably are initiated by larger M 23 C 6 carbides which are mainly precipitated on the grain and lath boundaries. The alignment of the dimple walls shows that fracture does not take place with void growth under Mode I, in which one would expect dimple wall elongation orthogonal to the surface, but a mixed Mode I and Mode II, where shearing is incorporated. The shearing direction acutely angled with the vertical is distinguishable from the fractograph of Figure 2c

Material Parameters for P91 Steel
The elastic bulk modulus of P91 at room temperature is K = 175,000 MPa and the shear modulus is G = 80,769 MPa. The isotropic hardening plasticity parameters are given in Table 1. Since no rate sensitivity is observed for P91 in the conditions of interest, the corresponding parameters are arranged to suppress the rate effect. Letting d x , d y and d z denote average inclusion diameters in the respective directions and S(%) and Mn (%) represent the weight percentages of sulphide and manganese in the matrix, respectively, the initial porosity can be estimated using quantitative metallography from Franklin's formula [25].
Assuming approximately spherical inclusions (d x = d y = d z ), the initial porosity is calculated as f 0 = 0.00044, based on the previously given chemical properties of P91. A higher value of f 0 was suggested ( f 0 = 0.002) in [26] thus in further sections the effect of f 0 on simulation results was investigated by using f 0 = 0.00044 and f 0 = 0.002. Franklin's formula serves as an estimate for the initial void volume fraction f 0 . The authors' simulations on unnotched disks show that initial porosity estimates within (0.0044, 0.002) resulted in only marginal differences in the computed failure times. Regarding the void nucleation, depending on the bimodal distribution of particles, it is possible to postulate two separate strain dependent probability distribution for each nucleation source. This requires mode detailed analysis which are not available to the authors. Hence, for the void nucleation kinetics we content with using previously suggested material parameters in the literature. In this context, we believe that separation nucleation mechanisms requires further studies. The parameters q 1 = 1.5, q 2 = 1, q 3 = q 2 1 of the extended Gurson model are chosen following [8,9]. Motivated by the fact that the volume fraction of the segregated inclusions f N is within a narrow band of 0.01 to 0.03 the parameters controlling void nucleation are chosen as f N = 0.02, N = 0.3 and S N = 0.1 [27,28]. The proposed range of k w for structural alloys is reported as 0 < k w < 3, see, e.g., [11]. In the current study, k w = 0 gives the best results for the predicted fracture strains. Nevertheless, a sensitivity analysis is performed investigating the effect of k w for various disk thicknesses.
Pertaining to the void coalescence and final fracture, the European Structural Integrity Society round robin [29] recommends the slope of the tail of the bilinear coalescence function (6) as 4. For the final void volume fraction at failure different references give different results, e.g., [30] takes f f = 0.25 in accordance with [31] whereas f f = 0.2 is used in [27,28,32,33], where the last two studies refer to steel A533 B. Based on the room temperature parameter fitting studies for P91 steels the coalescence and failure porosity are taken as f c = 0.1 and f f = 0.25.
Applying a nonlocal regularisation to local formulation of porous plasticity, introduction of a characteristic length scale is necessary. This length scale has been related to a physical quantity such as (four times) void size or (half) void spacing for ductile fracture mechanism [18]. If two populations of second phase particles are present, which is the case for P91, one should select the population which is dominant in crack initiation and propagation. While Xia and Shih observed that large inclusions constitute the main contribution to void formation while small inclusions only assist the hole link-up [34]. Thus, they selected the length scale as the mean spacing between voids nucleated by large inclusions. As seen from the fracture surface in Figure 2c, the small voids nucleated by MX precipitates are dominating in quantity and the mean distance between them ranges from 1.5 µm to 2.5 µm. On the other hand, the distance between large voids nucleated by M 23 C 6 precipitates are between 20 µm and 30 µm. Regarding the two studies mentioned, we use R = 5 µm as an average value considering the two population of precipitates. This value is also validated by the simulations.

Finite Element Model
2D axisymmetric simulations with CAX4R reduced integration elements for room temperature are conducted in ABAQUS/EXPLICIT with double precision. A solution of quasi-static problems with a dynamic-explicit solution procedure generally involves a very large number of time steps. In order to reduce the computational cost, mass scaling is applied with a target time step of 10 −3 s over the whole analysis which lasts for 150 s. This supplies acceleration of the simulations without changing the actual time scale of the process. The dies and the punch are modeled as rigid bodies and the disk as a deformable body. The interaction in-between the rigid and deformable bodies is assumed to be constant with a dynamic friction coefficient chosen as µ = 0.25, which is taken to be temperature independent. The friction coefficient has a direct influence on the necking as well as damage initiation location. In order to show this more clearly we have conducted a new set of simulations with and without friction effects. We have observed that the location of the localization gets closer to the center as the friction coefficient decreases. The thinning at the localization is more apparent with higher coefficients whereas it is less at the center. It is observed from the force-displacement curves that the crack initiation and failure takes place at lower displacements if the friction effect is higher due to the more pronounced thinning at the necking location.
The mesh used for 0.5 mm thickness is given in Figure 3. If using an integral-type nonlocal formulation, the element size should be below the characteristic length. For R = 5 µm which is selected depending on the microstructure of P91, this requirement puts severe restrictions on the selected mesh size. In order to reduce the computational time non-uniform meshing was used for the model: refined mesh at the crack location, coarse mesh at the clamped part. Initially, a uniform coarse mesh is applied to identify the approximate crack location. Then, refined meshing is applied to this region of interest.  Crack propagation is modeled using an element erosion technique: elements with Gauss points whose damage reaches the corresponding failure thresholds are removed from the computational stack. Although the applied nonlocal formulation removes the mesh dependence of the field distributions to an extent, for the crack propagation this could be limited. Hence, since element erosion is used, the direction of the crack and the crack propagation size are controlled by the mesh. Once biased and/or coarse mesh is used, predicted crack paths will be inevitably poor. In this study, sufficiently small sized elements with irregular distributions are used. The irregularity of the mesh created by the advancing front quad method accounts for the microstructural heterogeneity. With the same motivation the number of elements within the effective radius of each element vary spatially.

The Influence of Geometrical Parameters: Thickness
During SP test as the puncher deforms the SP specimen, multi-axial stress and strain fields occur. Various deformation stages take place: elastic bending, plastic bending, membrane stretching and, eventually, crack initiation and failure of the specimen.
The von Mises stress plots (Figure 4) of the SP disk with 0.5 mm thickness reveals that with the onset of contact between the puncher and the disk, high stresses occur underneath the contact region and this zone expands with the contact region till the crack initiation. Due to high stresses developed along with the contact, plastic deformation takes place underneath the contact region at the beginning of the test. As the puncher continues to deform the disk, maximum plastic deformation moves to the bottom surface and localizes in this region where crack initiates.  The deformation modes are intimately dependent on the specimen thickness. For thinner specimens, considerable thickness change and stretching occurs during the test until fracture. On the other hand, indentation caused by the localized plastic yielding underneath the puncher gets more pronounced with increasing thickness. This is revealed by the plastic strain plots of 0.2 mm, 0.5 mm and 1 mm disks at the initial stage (at displacement ∼0.07 mm and 0.3 mm) of the test given in Figure 5. The plastically deformed region under the contact area is more pronounced for the 1 mm thick disk compared to the thinner disks of 0.2 mm and 0.5 mm. As the puncher further penetrates through the disk, the plastically strained region for the 1 mm disk expands in the upper surface and moves away from the center, whereas for 0.2 mm and 0.5 mm it localizes at the bottom surface. In order to make further evaluations to reveal the influence of the disk thickness, we present the von Mises stress, plastic strain, pressure plots of 0.2 mm, 0.5 mm and 1 mm thick disks prior to crack initiation in Figure 6. For all thicknesses, maximum plastic strain occurs at the bottom surface of the disk where tensile stresses prevail and promote voidage. For the 1 mm disk there is also a highly strained region at the upper side close to the clamped region which occurs under the combined influence of tension and shear. Again for all thicknesses, the crack initiation starts from the bottom surface where plastic strain is maximum. Practically, the whole unclamped part of the disk exhibits the high stresses with a slight decrease at the crack initiation location for all thicknesses. If the hydrostatic stresses are compared (which is the negative of the pressure plotted in Figure 6c, for all thicknesses the unclamped region is under positive hydrostatic stress except the zone underneath puncher for 1 mm. As the thickness gets smaller, the distribution of the hydrostatic stress becomes more uniform through the thickness. Obviously, this has direct consequences on the void nucleation and growth. Figure 7 depicts the damage fields: f , f n , f g . Since in this figure k w is taken as 0, there is no damage growth due to shear. Hence, f is equal to the summation of f 0 , f n and f g . If pressure plots in Figure 6c are taken into account, one can see that just before crack initiation, nucleation mostly reaches its maximum value of 0.02 at the region under positive hydrostatic stress. Hence, with plastification all the void nucleation sources are exhausted where the distribution is rather uniform over the positive hydrostatic stress region. With thinner specimens the source of positive hydrostatic stresses is the stretching behavior. Thus underneath the punch there is considerable nucleation. This is not the case for increased thicknesses, where under the punch since the indentation mode is dominant, high compressive stresses with this confinement allows neither nucleation nor positive growth. f g normal also occurs in the region where positive hydrostatic stress prevails and has its maximum value where plastic strain is also maximal. For the disk with 1 mm thickness, the zone underneath the contact region where hydrostatic stress is negative f g is also negative. Thus, in this zone void shrinkage takes place instead of void growth. For all thicknesses, the crack initiates from the highest strained location underneath the puncher at the bottom surface and propagates upwards as shown in Figure 8, where plots of plastic strain are given in various stages of crack propagation for disks with thickness of 0.2 mm, 0.5 mm and 1 mm. This reveals the ductile crack propagation in the direction of the maximum plastic strain. Virtually maximum plastic strain occurs in the path due to localization of voidage and the crack simply follows this path since the main controlling variable for the crack path is the porosity. For the 0.2 mm disk, the final crack has an inclined path whereas for 0.5 mm disk the crack also starts with an inclined path but then kinks away almost half way. Similar kinking is also observed for this disk in the optical image of the section (see Figure 13c). As to the 1 mm disk, the crack initiating at the bottom side propagates vertically inside the shear band. Although there is a highly strained region at the upper side of the disk, both the maximum strain and damage localization occur at the bottom side where the crack initiates. In general, for thinner specimens the crack propagation direction is less inclined whereas for the thick ones, vertical crack propagation takes place. It is also seen that for the thin specimens relatively uniform fields are observed which is not the case for the thick ones.

The Influence of Material Parameters
The effect of shear softening due to void distortion and inter-void linking is investigated by using different values of the scaling parameter k w : k w = 0 and k w = 1 respectively representing no and full shear influence on ductile damage. These values are investigated for the two extreme thickness values: 0.2 mm and 1 mm. As concluded from Figure 9 increasing k w values results in an earlier loss of load carrying capacity, and thus, a decrease in the recorded displacements at incipient fracture. As anticipated, in the thinner disks the membrane stretching mode of deformation governs while shear effects are much less important. Accordingly, with an increase of thickness shear effects govern and as a consequence, the shear extension in the GTN model results in a considerable reduction in the materials ductility.
In Figure 10 the damage distribution plot at a displacement of 0.7 mm reveals how the k w parameter influences the total damage development in a 1 mm thick disk. In agreement with abovementioned observations, with the increase of k w from 0 to 1, the maximum damage value increases by 28% at the same displacement.
The effect of f 0 on force-displacement curves is investigated by taking f 0 as 0.002 from [29] and 0.00044 calculated from Franklin's formula. As anticipated, increased value of f 0 results in earlier fracture, while the hardening part of the curve is almost not affected (see Figure 11a). In Figure 11b,c the damage development of a 0.5 mm thick disk prior to crack initiation for f 0 = 0.002 and f 0 = 0.00044 are presented showing that with higher initial porosity the damage accumulation also increases.
min max  In conclusion, if the shear extension is utilized, a premature fracture occurred, so the k w parameter had to be taken as 0 to achieve agreement with experimental results. As to the initial porosity, while the effect on force-displacement curves was not so prominent, the f 0 = 0.00044 as calculated from Franklin's formula gave closer agreement with the experimental curve. Additionally, SP tests for largest and smallest thicknesses with and without damage effects were simulated to see its effect on necking. In all cases it is observed that necking occurs with or without damage whereas with damage its evolution with thinning until the final rupture occurs much faster.

Comparison to Experimentally Determined Results
In Figure 12 the computational and experimental load-displacement curves of SP disks are compared. Experimentally two to three tests were carried out for each thickness. For comparison with the numerical curves, averages of experimental load-displacement curves are used. The results show a good agreement between computational and experimental curves especially in terms of hardening. The maximum strengths are slightly underestimated for disks with 0.7 mm, 0.8 mm and 1 mm thickness. Generally speaking, the calculated curves exhibited a steeper force drop compared to the experiments due to axisymmetry assumption in the simulations which prevails in both localization and fracture behaviour. It should be noted that for the selected sample size, shape, loading conditions and selected mesh size/length scale ratio both local and nonlocal formulation estimations were similar. The similarity was revealed in terms of force displacement plots, i.e., energy dissipation during fracture, fracture patterns and was attributed to the milder stress and strain gradients. Considering the severe change in element aspect ratios under stretching during SP test simulation, only nonlocal formulations supply spatially invariant material length. With this property, the developed framework constitutes a unified modeling environment for problems involving sharp notches where high stress gradients are evident or problems with smaller specimen sizes where size effects will govern.  Table 1. Discrete circles represent the experimental data and the solid lines represent the simulations.
In Figure 13 a comparison of the deformed profiles obtained from optical microscopy observations (left) and simulations (right) is given. The fractured SP disks are sectioned by a precision cutter and images are obtained by an optical microscope. Both crack locations and the paths of the cracks of the computational results are found to be in good agreement with the experimental ones. It is noted both from FE results and optical images that the location of crack initiation slightly moves further away from the center of the disk with increasing thickness. In 3D, the crack follows approximately a circular path along the maximum plastic strain contour for all thicknesses (see Figure 14) which complies with the axisymmetry assumption applied in the present modeling.

Conclusions and Outlook
In this study, a combined experimental and computational investigation of deformation and fracture during room temperature SP tests of P91 steel disks has been presented. To this end, tests were conducted on disks with different thicknesses. A void driven ductile fracture mode is recorded for all cases which is preceded by a diffuse necking with membrane stretching followed by a localized deformation for smaller thicknesses and shear localization for larger ones.
On the computational side a rate dependent porous-plastic constitutive model with a non-local extension is established to predict the deformation and fracture behavior of P91 steel during SP testing. The nonlocal formulation allows a natural control of the localization size through incorporating a material length parameter relating to the inter-particle spacing. With the developed framework besides a detailed full field investigation of the deformation process in the SP test, parameter sensitivity analyses are also realized. The set of material parameters were identified through a combined quantitative metallurgical and inverse mechanical analysis by comparison of the numerical and experimental force-displacement curves. The accuracy of the simulation results with the identified parameters were assessed by comparing experimental and computational force-displacement curves as well as the deformed sections of disks obtained from FE and optical images. A good agreement is captured between the experimentally and computationally determined force-displacement curves as well as the deformed sections at fracture.
The parameter sensitivity analysis conducted for the shear damage parameter k w shows that increasing values of k w result in earlier fracture and this effect got more pronounced with increasing disk thickness due to the influence of shear. To a lesser extent initial porosity, f 0 , also affected the force-displacement curves in the same way since with higher values of f 0 the damage development was higher. The developed framework allows modeling ductile fracture in the presence of sharp stress gradients driven by, e.g., notches. The model was also found to be efficient in predicting the deformed geometry: necking patterns, crack initiation location and the crack propagation paths by comparing the sections of the failed samples obtained through simulations and optical microscopy.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: SP small punch GTN Gurson-Tvergaard-Needleman RKR Ritchie-Rice-Knott YS yield stress UTS ultimate tensile strength SEM scanning electrone microscope