A Compensation Algorithm for Large Element Characterizing the Damage Evolution Process and Its Application to Structure Collisions

When simulating the process from elastic–plastic deformation, damage to failure in a metal structure collision, it is necessary to use the large shell element due to the calculation efficiency, but this would affect the accuracy of damage evolution simulation. The compensation algorithm adjusting failure strain according to element size is usually used in the damage model to deal with the problem. In this paper, a new nonlinear compensation algorithm between failure strain and element size was proposed, which was incorporated in the damage model GISSMO (Generalized incremental stress state dependent damage model) to characterize ductile fracture. And associated material parameters were calibrated based on tensile experiments of aluminum alloy specimens with notches. Simulation and experimental results show that the new compensation algorithm significantly reduces the dependence of element size compared with the constant failure strain model and the damage model with the linear compensation algorithm. During the axial splitting process of a circular tubular structure, the new compensation algorithm keeps the failure prediction errors low over the stress states ranging from shear to biaxial tension, and achieves the objective prediction of the damage evolution process. This study demonstrates how the compensation algorithm resolves the contradiction between large element size and fracture prediction accuracy, and this facilitates the use of the damage model in ductile fracture prediction for engineering structures.


Introduction
In a collision, engineering metal structures deform severely and fracture may occur locally. Here, fracture occurs under a monotonic load, which is different with fatigue fracture that damage is accumulated under repeated loads [1]. And fracture is determined by material properties and the stress state. Mechanical properties vary with materials, e.g., welding parameters would affect the tensile strength [2]. For the stress state, shear fracture happens with low stress triaxiality and necking fracture with high stress triaxiality [3]. To simulate the phenomenon, three material constitutive models are required: an elastic-plastic model characterizing the elastic-plastic deformation; a damage initiation model characterizing the occurring of damage; a damage evolution model characterizing the material degeneration due to damage. However, element size is not consistent for engineering structure and constitutive models.
Due to the limitation of computational capability, large shell elements are usually used for deformation of engineering structures. While small solid elements are needed for constitutive models, especially for damage models. This is because fracture generally happens in the small region, thus the deformation gradient is very large, which demands small solid elements to describe the severe deformation accurately [4,5]. Therefore, for a fracture prediction problem of engineering structure with large shell elements, there are two key issues to be evaluated: the stress state absence of shell elements compared with solid elements and the effect of element size on constitutive model. Among them, the relationship between element size and the damage model is studied in this paper.
In the field of constitutive model, Swift and Voce models are popular isotropic hardening rules, reflecting the relationship between stress and plastic strain [5,6]. For the damage initiation model, the forming limit diagram [7] (FLD) can be used in proportional loading situations; in non-proportional loading situations, the FLD can be converted to the relationship between critical strain and stress state, based on which the nonlinear accumulation algorithm is applied [8]. Since the FLC mainly focuses on necking failure and does not pay enough attention to shear failure [9], the modified Mohr-Coulomb criterion (MMC) [10] characterizing the relationship between critical strain and stress state can be used. In the principal strain space, compared with the FLC curve, there is an additional critical curve in the MMC criterion that considers the stress state from uniaxial compression to uniaxial tension, which includes shear failure [5].
Damage evolution models are divided into two types: meso and macro. The mesoscopic damage model focuses on the evolution of mesoscopic material defects, and the most representative model is the Gurson-Tvergaard-Needleman model (GTN) [11]. And the macroscopic damage model directly characterizes the material load-bearing degradation phenomenon through damage variables, the most representative model is the coupling relationship between damage and stress proposed by Lemaitre [12]. Neukamm et al. proposed a generalized incremental stress-state dependent damage model (GISSMO) in 2008, using the Lemaitre-type coupling of damage and stress to realize the damage evolution process under complex stress states [8,13,14]. Then the model was widely used to study the failure of metal structures in forming, compression and joining situations [15][16][17].
In the field of test design, ductile fracture under various stress states is achieved by changing the loading state of the failure region. For bars, the fracture of stress triaxiality ranging from − 0.33 to 0 can be achieved with specimens of different height-to-diameter ratios under compression condition [18]. The fracture of stress triaxiality ranging from 0.33 to 1 can be achieved with arc-notched specimens under tensile condition, and the stress triaxiality of the failure area can be adjusted by changing the arc notch radius of specimens [19,20]. For sheets, which are prone to transverse instability under compression condition, generally only tensile tests are conducted to achieve fracture of stress triaxiality ranging from 0 to 0.66 by changing the angle between the minimum section of specimens and tensile force. Regarding the realization of the angle, one way is to design a series of specimens with different notches [21]. Because fracture usually originates from notches under loads, desired stress states can be obtained from specimens with notches [22]. And another way is to keep the specimen unchanged and change the loading angle by means of clamps [23,24].
In the field of ductile failure simulation with shell elements for engineering structures, owing to the serious deformation close to the fracture region, measure strain increases nonlinearly with gauge length decreasing in the experiment, and the relation between calculated stain and element size in the finite element (FE) model is also similar to this. To ensure the objective of the FE solution, a compensation term related to the ratio of thickness to size of the shell element is usually adopted in the damage model. Alsos et al. used the compensation term in the RTCL damage model and simulated the penetration of stiffened panels, which matched well with experiments [25,26]. Storheim et al. also use the compensation term in the damage model [27]. Moreover, some studies found the effect of element size on damage accumulation is also relevant to the stress state. Walters proposed a relation among failure strain, the ratio of thickness to size of shell elements and the stress triaxiality for the damage model [28]. Korgesaa et al. used the relation in the ductile fracture simulation of stiffened panels [29]. Korgesaar also used the relation in the RTCL damage model to simulate the collision failure of the hull at low stress triaxiality [30]. And some conclusions were found: The effect of element size in shear state is less than that in tensile state [31]. Uniaxial tensile state is more sensitive to damage accumulation with respect to element size than plane strain and biaxial tensile states [32].
In the impact accident, the kinematic energy is usually dissipated by the deformation of energy-absorbing structures. They are generally thin-wall structures, including circle tubes [33], square tubes [34], honeycombs [35], etc. In recent years, the splitting of the tube with notches, owing to the large stroke to length ratio and high specific energy absorption, is employed in the energy-absorbing structures [36]. And the simulation method is usually used to obtain the optimum structure parameters [37], thus the simulation accuracy is very important. But the splitting process is difficult to simulate accurately. This is because the splitting process involves deformation concentration caused by notches, material degeneration and crack growth due to void growth, the effect of stress states, etc. Guan et al. used a Johnson-Cook failure model to predict the splitting of circle tubes, where fracture variation dependent on the stress state was considered, but the material degeneration derived by void growth was ignored [37]. Rouzegar et al. used a surfacebased cohesive technique to model the splitting of circle tubes, where material degeneration was considered, but the fracture paths were determined before simulation [38]. Owing to the complexity of simulating the splitting of a circle tube, this was chosen to check an algorithm relevant to damage evolution proposed in this paper.
In this research, a compensation algorithm for large element charactering the damage evolution process was proposed, which establishes a nonlinear relationship between element size and failure strain, and the ductile failure simulation of thin-walled metal structures under collision loading was performed. First, specimens of aluminum alloy 5083-O sheet with different notches were designed to achieve fractures of different stress states. Secondly, the Lemaitre-type coupling relationship between stress and damage was established for the gradual deterioration of bearing capacity exhibited during material ductile fracture. Then, the compensation algorithm for large element was incorporated, and the calibration of the model parameters was completed by a surrogate model and an intelligent optimization algorithm, and objective simulation of the experiments was achieved. Finally, the damage evolution simulation of axial splitting of a circular tube was realized.

Experiments and Analysis of Aluminum Alloy 5083-O Notched Specimens
In order to achieve ductile fracture due to different mechanisms under plane stress state, the standard tensile specimen, the specimen with arc-shaped notches and the specimen with sharp-angled notches shown in Figure 1 were designed. The detailed dimensions are shown in Table 1, the thickness of the notched specimens is 2.5 mm, and the thickness of the standard specimens is 2 mm. The thickness difference of specimens was caused by adjustment during processing. The specimens were cut by wire EDM along the rolling direction on an aluminum alloy 5083-O sheet of 1000 × 1000 mm 2 . The "degree" in the specimens with sharp-angled notches refers to the angle between the tangential direction of the smallest section in the specimen plane and the tensile direction, expecting that as the angle increases, the ductile fracture mechanism changes from shear fracture to necking fracture. Under tensile loading, the deformation of specimens is concentrated around the smallest section, and specimens are not subjected to force in the out-ofplane direction, which is plane stress state. Tensile experiments were performed at room temperature using an INSTRON 8801 fatigue test machine, and the displacement of specimens was measured using a video extensometer with the gauge length 50 mm. The tensile speed was 3 mm/min. To maintain the accuracy of experiments, the average of three repeated experiments was used.
The force-displacement curves of experiments are shown in Figure 2. All curves have three significant phases, the straight line segment before the circle marker is the linear elastic phase, and the curve after the marker is the plastic deformation phase and the damage evolution phase coupled with it. The damage initiation point needs to be determined by calibrating the constitutive model. And the rapidly falling part indicates the fracture response, which corresponds to the failure strain in the constitutive relationship. The peak force and failure displacement of different types of specimens are significantly varied, which is mainly caused by the excessive difference in geometric configuration and the different area involved in deformation. The deformation response difference of the same type of specimens is mainly from the effect of stress state on the damage evolution. In the response of the specimens with arc notches shown in Figure 2b, as the notch radius increases, the specimen stress state is closer to uniaxial tension and the failure strain gradually increases. In the response of the specimens with sharp-angled notches shown in Figure 2c, as the angle between the tangential direction of the smallest section in the specimen plane and the tensile direction increases, the specimen stress state changes from shear to tension and the failure strain gradually decreases.

Constitutive Model Based on Strain Element Homogenization and Calibration
Corresponding to the experiment curves, the elasticplastic model in this paper is as follows: where σ and σ s are the equivalent stress and the yield stress respectively, ε and ε p are the equivalent strain and the equivalent plastic strain respectively, E is the elastic modulus, k 0 , Q and γ are constants from the Voce hardening law [5], and k 0 = σ s . To account for non-proportional loading paths, in the damage initiation model the instability variable F is defined as [14]: where ε c (η) is the critical equivalent strain related to the stress triaxiality η, σ m is the mean stress.
The relationship between the critical equivalent strain and the stress triaxiality represented by the MMC criterion is defined as follows [5,10]: where A and n are the material constants related to the plastic deformation, and c 1c , c 2c , c 3c are the material constants for the damage initiation.
When the instability variable F in Eq. (2) reaches unity, the damage is coupled to the stress according to the GISSMO model and the coupling relationship is as follows [14]: where σ * and σ are the equivalent stress with and without damage respectively, D and D c are the current damage value and the damage value at the onset of stress-damage coupling respectively, m is the material constant for the effect of damage on stress. When plastic deformation occurs, the damage starts to accumulate by Eq. (9): , where ε f (η, l e , t) is the equivalent failure strain which is related to element size l e , thickness t and stress triaxiality η. When damage value D reaches unity, the material fails. Based on the relationship between the equivalent failure strain and element size from Ref. [28], this paper introduces a calibration parameter regarding the reference element size by Eq. (10): where ε c (η) is the critical equivalent strain in Eq. (4), l e , and t are the element size and the thickness of shell elements respectively, S is the area of shell elements, ε f0 (η) is the equivalent failure strain of the reference element related to triaxiality η, l e0 and t 0 are the element size and the thickness of the reference shell element respectively, r is the parameter related to element size.
It's worthwhile to mention that no matter what shape, e.g., square, rectangle and triangle, the element is, l e represents the equivalent size as the area keeps constant.
The relationship between the equivalent failure strain and the stress triaxiality for the reference element is still characterized by the MMC criterion Eq. (4). The values of parameters c 1c , c 2c and c 3c are different from the damage initiation model and are denoted by c 1f , c 2f and c 3f now, while the values of A and n are consistent with the damage initiation model.
The calibration of the model parameters was done using a surrogate model and an intelligent optimization algorithm. The calibration flow is shown in Figure 3. Sample points are determined based on the D-optimality design criterion. And the surrogate model is a linear polynomial response surface. The optimization based on the surrogate model includes two steps: first, use the adaptive simulated annealing algorithm to find an approximate global optimum; second, starting from the solution, use the gradient-based dynamic leap-frog method to find the local optimum. Thus, the optimum is obtained. Further, the domain of parameters is adapted based on the accuracy of the previous optimum [39]. The calibration process was done in the LS-OPT.
As shown in Figure 3, the calibration includes two parts: the elastic-plastic model and the failure model. As we all know, the elastic-plastic model has no relation with specimen geometry, thus all specimens were adopted to calibrate its parameters. However, the failure model is related to element size. A general method is calibrating failure model with reference element firstly, then calibrating element size regularization parameters with other elements [14]. In this research, the shell element is used, thus two parameters, thickness and element size, need to be concerned. The element with thickness 2.5 mm and size 0.2 mm is employed as the reference element, and specimens with thickness 2.5 mm, namely specimens with notches, are adopted to calibrate failure model parameters m, A, n, c 1c , c 2c , c 3c , c 1f , c 2f , c 3f . Then, element size 0.3 mm, 0.5 mm are used to calibrate element size regularization parameter r. Since the regularization function (10) takes thickness and element size into account, the standard specimen with thickness 2 mm and notched specimens with thickness 2.5 mm are both adopted in the calibration of parameter r. As to the determination of element size, both efficiency and accuracy are considered, which is described with simulation results in Section 4.

Simulation and Discussion
The FE models of the specimens in Figure 1 were established using fully integrated shell elements with three integration points along the thickness direction. In order to compare the effect of the constitutive model on failure prediction, the constant failure strain model (CFS), the damage model based on the linear relationship of Ref.
[28] (r = 1) and the damage model based on the nonlinear relationship Eq. (10) (r = 2.39) were used. The material The elastic-plastic model      Tables 2  and 3. Among them, the constant failure strain is the elongation at break of the ST specimen. Three element sizes of 0.2 mm, 0.3 mm, and 0.5 mm were selected and 0.2 mm was used as the reference size, which was determined considering the accuracy and efficiency of calculation. With the strain gradient increasing, element size should be decreased to describe the deformation accurately. Among specimens, the maximum strain gradient appeared in the specimens with sharpangled notches, and they were used to determine the appropriate element size. The N-45 specimen was chosen to analyze. The force-displacement curves are shown in Figure 4c, we can see that the curves of element size 0.2 mm and 0.3 mm almost coincide before degeneration. That means element size 0.2 mm is enough to describe elastic-plastic deformation of specimens. And if employ smaller elements, it is believed that the accuracy won't increase significantly while the computation efficiency will decrease a lot. So 0.2 mm was chosen as the reference size. It's worthwhile to mention that damage behavior is not considered in the determination of the reference element size. This is because the damage model is relevant to element size. If the element size is appropriate according to the plastic model, the damage model should have the capability to characterize the material degeneration accurately. Nine simulation results were obtained for each specimen model. And the force-displacement curves comparison of ST specimen, N-R10 specimen, and N-45 specimen in different simulation condition is shown in Figure 4. Obviously, the simulated force-displacement curves of ST specimen and N-R10 specimen are in good agreement with the experiments, and the failure model has less influence on the simulation results, and the results of ST specimen are better between them; while the simulated force-displacement curves of N-45 specimen in Figure 4c are highly correlated with the failure model.
The fracture specimens and the equivalent plastic strain nephograms at fracture using the nonlinear relationship Eq. (10) with element size 0.2 is shown in Figure 5. Fracture locations are similar for both experiments and simulations. The plastic deformation of ST specimen in Figure 5a is uniformly distributed in the parallel region while the damage failure occurs locally with a very small percentage in the plastic deformation region. The force-displacement curve in Figure 4 characterizes the overall deformation of the specimen. For ST specimen in Figure 4a, the force-displacement curve mainly reflects the plastic deformation and cannot reflect the difference between different failure models. N-R10 specimen in Figure 5b is similar with ST specimen, except that the plastic deformation region of the specimen is relatively small, thus the damage failure region has a relatively high proportion in the plastic deformation region, so the force-displacement curves in Figure 4b  the difference of different failure models to a certain extent. The N-45 specimen in Figure 5c is exactly the opposite of them, and the damage failure region basically overlaps with the plastic deformation region, thus the force-displacement curves in Figure 4c reflect the difference between different failure models clearly. From the force-displacement curves of N-45 specimen in Figure 4c, it can be seen that the 3 solid lines as to r = 2.39 of different element sizes coincide well and are in good agreement with the experiment, the 3 long clash lines as to r = 1 of different element sizes coincide poorly, and the 3 short clash lines as to CFS model of different element sizes coincide worst. Due to strain element homogenization, the simulation results of CFS model relate to element size. The correlation of the simulation result with element size is significantly reduced after using compensation algorithms, and the nonlinear algorithm (r = 2.39) works best between them, which effectively guarantees the objectivity of the failure prediction.
The formula for calculating the simulation error is where F exp and F sim are the experiment force and simulation force corresponding to the same displacement s respectively, s max is the maximum displacement observed in the test and simulation. In addition, all elements on the fracture surface was used to calculate the stress triaxiality of specimens. Because the fracture volume between simulation and experiment is different, the geometry of the FE model differs from that of the specimen once an element is deleted, and the average triaxiality before element deletion is adopted to represent the stress state of the specimen. The stress triaxiality values of specimens calculated based on the results of element size 0.2 mm and the failure model r = 2.39 are listed in Table 4, and the stress triaxiality values based on other constitutive models and other elements are similar to these. Thus the simulation error calculated based on 5 specimens with sharp-angled notches and different failure models are shown in Figure 6, obviously the error of r = 2.39 stays low over the stress states ranging from shear ( η = 0 ) to biaxial tension ( η = 0.66 ), while the error of the CFS model is high overall, and the error of failure model r = 1 varies with the simulation condition. When simulating the shear stress state with large shell elements (0.5 mm), the error of r = 1 is the largest, approaching 100%; but when the stress state is close to the biaxial tension, the error of r = 1 is the smallest, almost 0. In addition, the error of the CFS model deceases with element size, this is because the fracture energy increases with element size, and when using the element sizes 0.2 mm, 0.3 mm, 0.5 mm in turn, the fracture energy increases and approaches that of experiment, which can be seen from Figure 4c. But with larger elements, fracture energy may excess that of experiment, and the error will increase with element size. Furthermore, the effect of element size on fracture simulation should also be analyzed deeply. The failure mechanism between simulation and experiment is different. In physics, fracture appears with cracks propagation, which occurs in a small region. While in FEM, fracture occurs with the successive deletion of elements, the volume of which is dependent with element size and is generally larger than that of cracks. Thus, the experiment force decreases smoothly, but for the CFS model response in Figure 4c, simulation force decreases with a stair-step shape consistent with element deletion, that is there is a force maintenance and even hardening and decreasing stage for each deleted element. And due to failure strain is constant, fracture energy is consistent with volume of deleted elements and increases with element size, that is the area below the force-displacement curve increases Table 4 Stress triaxiality values of specimens with element size. To ensure the objective of simulation, a compensation algorithm for failure strain is adopted, e.g., Eq. (10), failure strain decreases with element size, which can also be seen in Figure 7, the equivalent plastic strain nephograms at fracture for N-45 specimen using the model r = 2.39 with several element sizes. Thus, we can see from Figure 4c the force discrepancies between simulation and experiment for the model r = 2.39 is much smaller than that of the CFS model. However, with element size increasing, the discrepancy of failure volume between simulation and experiment still increases. The compensation algorithm is just a modification for failure strain, which can only reduce the discrepancy of overall response to some extent.
To guarantee the objectivity of the damage model, the nonlinear relationship between failure strain and element size shown in Eq. (10) was proposed. The variation of failure strain with element size at a certain stress state (stress triaxiality 0.33) is shown in Figure 8, and it can be seen that: (1) The failure strain gradually increases with the increase of t/l e and gradually tends to a constant value with the decrease of t/l e .
(2) The parameter r characterizes the nonlinear relationship between the failure strain ε f and t/l e . When r = 1, the relationship is linear, which is the same with Ref. [28].
(3) The intersection point A of all curves characterizes the failure strain corresponding to the reference element.
From the simulation error in Figure 6, it can be seen that the correct calibration of parameter r based on material property can effectively improve the simulation accuracy and avoid the less accurate prediction of linear model r = 1 for certain simulation conditions.

Damage Evolution Simulation of Axially
Splitting of a Circular Tubular Structure

FE Model
The combined splitting circular tube energy absorber is used to dissipate collision energy of rail vehicles [37,38], and the structure sketch is shown in Figure 9, consisting of a circular tube and a die. The thickness of the circular tube is 2 mm, and 6 notches are evenly distributed at the bottom of the circular tube. These notches are used to guide fracture paths due to the weakness around notches [22]. The friction coefficient between the tube and the die is 0.2. To model impact process, a quasi-static downward compression of 40 mm is applied to the top of the circular tube. And lateral displacement constraints are also applied to the top of the circular tube to ensure the desired deformation mode. For the die, rigid material is adopted, and all displacements are constrained. The FE models were established using fully integrated shell elements with three integration points along the thickness direction. 4 element sizes 0.2 mm, 0.5 mm, 1 mm, 2 mm and 3 failure models the constant failure strain model (CFS), the damage model based on the linear relationship of Ref. [28] (r = 1) and the damage model based on the nonlinear relationship Eq. (10) (r = 2.39) were used for the tube. The compression simulations were completed in the LS-DYNA with an implicit solver. It's worth mentioning that these element sizes were determined according to the elastic-plastic deformation of the structure. The force-displacement curves of r = 2.39, r = 1, CFS models are shown in Figure 10, the curves before first peak represent the elastic-plastic deformation, obviously, the deformation of element size 0.2 mm matches well with that of element size 0.5 mm, and the deformation of element size 1 mm matches well with that of element size 2 mm. The deformation difference between element size 0.2 mm and 1 mm is caused by the simplification on the geometry for large element, which will be discussed later. And the first peak force difference is caused by failure models. So, these element sizes are appropriate to characterize the elastic-plastic deformation of the structure. Moreover, damage behavior is not considered in the determination of element size, which is the same with that of the reference element size discussed in Section 4.
Further, t/l e for the structure is 10, 4, 2, 1, some of which are below the ratio range of the model calibration (12.5, 8.33, 5). From ratios, with the ratio decreasing, failure volume in simulation increases, and the discrepancy between simulation and physics increases, thus the simulation accuracy would decrease. Therefore, the simulation accuracy of the structure would be smaller than that of specimens.

Deformation Analysis and Comparison of Damage Models
The deformation modes under different conditions are similar, and the result of element size 0.2 mm and failure model r = 2.39 is analyzed as an example. The equivalent stress nephograms are shown in Figure 11, including 4 deformation stages: (1) In the initial stage of the compression, the circular tube starts to contact with the die and the impact force rises sharply, reaching the peak impact force of 5.23 kN at the displacement of 3.2 mm.
(2) Then the tube wall bends, the material around the notches fails and the impact force gradually decreases, reaching the minimum value of 3.84 kN at the displacement of 7.6 mm.
(3) As the contact area between the circular tube and the die increases, the total friction force gradually increases and the impact force rises again, reaching 4.82 kN at the displacement of 11.6 mm.
(4) Finally, the structural deformation stabilizes and the impact force keeps constant.
The force-displacement curves of r = 2.39, r = 1, CFS model are shown in Figure 10, and there are some differences in the results of different element sizes, among which the CFS model is the worst case. Taking the results of the element size 0.2 mm as the benchmark, the errors caused by the element size calculated with Eq. (12) are shown in Figure 12, and the errors of CFS model is much larger than damage models. Apparently, the homogenization effect of the element on the failure strain has a significant effect on the force response, and the introduction of the compensation algorithm based on the damage model effectively reduces the correlation of the simulation results with the element size.
The energy dissipation mechanism should be studied deeply. The dissipated energy is derived from two parts: the plastic deformation of the tube and the fracture energy of the material around notches. And the plastic deformation dominates because the energy dissipated by the failed elements only accounts for 17.2% of the total internal energy according to the simulation result of the nonlinear algorithm r = 2.39 using the elements of size 0.2 mm. The extent of plastic deformation can be characterized by the average equivalent strain calculated by the equivalent strains at the 6 middle points. The middle points are determined according notches, e.g., in Figure 9, the point B is a middle point for the undeformed structure. The equivalent strain is relevant with the equivalent stress, which is determined by the fracture tolerance around the notches. That is if the material around the notches is easy to fracture, the equivalent stress and the equivalent strain at the middle point between two notches are small, then the compressive force is small and associated energy dissipated by the structure is also small, and vice versa.
It's worthwhile to mention the steady force difference among different element sizes. The force at the displacement of 30 mm is used to study the steady force  difference, and the CFS model is taken as an example. From the equivalent plastic strain nephograms at the displacement 30 mm shown in Table 5, with element size increasing, the deformation concentration around these notches decreases, that means it is hard to fracture, thus the relevant average equivalent strains at middle points increases shown in Table 5, so the steady force increases with element size shown in Figure 10c. In addition, the stress-plastic strain curves of failure elements are shown in Figure 13, and the failure strain decreases with element size in the compensation algorithms, which means fracture is easy to happen as the element size increases. This reduces the effect of deformation concentration to a certain extent. Therefore, for the models r = 1 and r = 2.39, the steady force discrepancies caused by element size shown in Figure 10a, b are smaller than that of the CFS model shown in Figure 10c.
From Figure 10, we can see that there are abnormal peak force values at the displacement 7.6 mm for element sizes 1 mm and 2 mm. Owing to the large element size, the notch width 0.2 mm is ignored in the models with element size 1 mm and 2 mm shown in Table 5. The CFS model still is taken to analyze, the plastic strain nephograms at the displacement 7.6 mm are shown in Table 5. For element sizes 1 mm and 2 mm, almost 4 elements around each notch meet the failure criterion, while only 1 element fails for element sizes 0.2 mm and 0.5 mm. Obviously, compared with element sizes 0.2 mm and 0.5 mm, the failure regions for element sizes 1 mm and 2 mm are much larger, and this induces abnormal peak force values in these two situations. Large elements are always used to achieve high efficiency, but the induced geometry simplification may cause abnormal response. This should be given special analysis. For this case, although the abnormal peak force appears when using element sizes 1 mm and 2 mm, the steady force is similar with element sizes 0.2 mm and 0.5 mm in the nonlinear algorithm r = 2.39. When analyzing dissipated energy, since Eq. (12) actually evaluates dissipated energy error, the errors caused by element sizes can be seen from Figure 12. The maximum value is under 30%, and the error for element size 1 mm is under 15%. Obviously, the dissipated energy calculated from element size 1 mm and 2 mm can be accepted. But if analyzing peak force, the results from large elements are apparently wrong. In a word, whether calculation from large elements can be accepted depends on what index to concern. From force-displacement curves shown in Figure 10 and errors shown in Figure 12, we can see that the difference between compensation algorithms r = 2.39 and r = 1 is small. Considering the fact that plastic deformation dominates the energy dissipation mechanism as discussed previously, the reducing effects of the two compensation algorithms on deformation concentration caused by element size are similar, which results in similar force responses. However, for failure elements, there is a big difference in the stress-strain curves between the two compensation algorithms shown in Figure 13. That is fracture energy is different for the two compensation algorithms. Therefore, for the structure where fracture energy dominates energy dissipation mechanism, compensation algorithms would affect a lot on the results. The tube with relatively high thickness is one case. When it is compressed, plastic hinges are the main deformation mode. Owing to the high thickness, the fracture may occur on the hinges [40]. And this will be study in the next research.

Conclusions
(1) Ductile fracture under various stress states was experimented and simulated based on specimens with notches. In the simulation, fracture volume is positively relevant to element size, and the result of the constant failure strain model is seriously related to element size. To ensure result objective, the compensation algorithm adjusting failure strain according to element size is incorporated in the damage model to keep fracture energy constant. And the results show that fracture simulation dependency on element size reduces a lot. (2) A compensation algorithm reflecting the nonlinear relationship between failure strain and element size was proposed. The new algorithm keeps the fail-ure prediction errors of different element sizes at a lower level in all stress states, while the errors of the existing linear compensation algorithm vary with the stress state. (3) The damage evolution process of axially splitting of a circular tubular structure was simulated. Owing to the deformation concentration around notches is relevant with element size, the results of the constant failure strain model are seriously dependent on element size. And because compensation algorithms have similar reducing effect to the deformation concentration, both damage models with the linear compensation algorithm and the nonlinear compensation algorithm can objectively predict the ductile fracture process. And the fracture energy values difference between the two compensation algorithms just make accuracy of the nonlinear compensation algorithm better a little. (4) The structure where fracture energy dominates energy dissipation mechanism, e.g., the compression of a tube with high thickness, will be studied in the future, where the new nonlinear compensation algorithm would play an important role.