Simulation of four-point beam bending test using the X-FEM method

Funding: This research was financed by the Lublin University of Technology Science Financing Subsidy FN16/ILT/2020. Abstract: This article presents the results of the computer simulations of a four-point bending test of a concrete single-edge notched beam. In this publication, the authors compared the X-FEM method of simulating the crack in the Abaqus FEA system. The paper also contains the results obtained with the Abaqus subroutine recently developed by the authors, used for defining the direction of the crack and the failure criterion. The publication explains the way of working of this algorithm. The described calculations show that computer analysis gives unrealistic results in terms of the destructive force. SEN-beam is an interesting laboratory test in which a lot of factors influence the results. It is especially important to study what phenomena occur in the final phase of the study, when the crack tip is near the opposite side.


Introduction
This article presents an attempt to simulate a four-point bending test on a concrete singleedge notched beam (SEN-beam) in the Abaqus FEA system using different material damage methods, especially the X-FEM method with its own damage criteria. The aim of this article is to verify how the X-FEM method behaves when simulating the SEN-beam, and check the Abaqus User Subroutines developed by the authors. The SEN-beam test is an interesting topic, but there are not many publications about it. It was first described by Carpinteri et al. [1]. The paper [2] describes different tests including the SEN-beam bending test. It uses the cohe-sive-zone damage model with different approaches to implement this method in numerical calculations. The publication described here, however, does not focus on concrete, but only on the cohesive-zone method for any brittle material. Another work [3] describes different fracture modelling approaches to gravity dam failure. Various models are described, including a concrete SEN-beam. This paper is a good data source. The paper gives, for example, crack paths for various numerical methods as well as for laboratory tests of the described beam. The paper [4] by Tabatabaee et al. includes, among others, a comparison of the SEN-beam test with other similar tests.

Description of the model
Dimensions and loads of the test are taken from Schalangen's paper [5] and are shown in Fig. 1. Lengths are indicated in mm. Fig. 1. The geometry of SEN-beam. Source: the authors' own study based on [5] The external dimensions of the beam are 440x100x100 mm. In the centre of the bottom edge, there is a 20mm deep notch, from which the crack should start. The simulations were performed in the Abaqus FEA system. The load was modelled as displacement at two points, which increased linearly to the value of 0.909 mm at the right point and to the value of 0.091 mm at the left point. The load was applied cyclically (the model was loaded and unloaded 5 times during calculations), because the task has a certain characteristic -the force in the supports increases proportionally to the given displacement, even after the failure of the element. During the tests, it turned out that the crack length increases even when the force stays the same or decreases. Therefore, it is not appropriate to inflict a constantly increasing load in such a task, because the crack cannot keep up with the load, hence the need for cyclic unloading.
The model mesh is shown in Fig. 2. The mesh size is about 15 mm, while in the area where the crack is predicted, the mesh has been compacted to 1.5 mm. The taken material is concrete C35/45 with the following parameters: • Young's modulus -E = 34 GPa; • Poisson's ratio -ν = 0.2; • Compressive strength -f ck,cube = 45 MPa; • Tensile strength -f ctm = 3.2 MPa; • Critical strain energy release rate -G Ic = 0.1 N/mm. These parameters were taken from the papers [6] and [7].

Explanation of the X-FEM method
X-FEM (extended Finite Element Method) is a method of simulating the fracture regardless of the finite element mesh. In each increment, the program leads the crack to the next element, modifying the displacements in the element nodes with an appropriate shape function. The advantage of this method, for example, over the methods which use smeared crack or element deletion is a short calculation time (the longest calculations used in this test lasted 15 minutes), high accuracy in predicting the crack path regardless of the finite element mesh. This method requires little data to calculate the simplest maximum principal stresses criterion. The disadvantage of this method is that it is limited to one crack. What is more, it works best only for static calculations.
The authors attempted to simulate the fracture of the SEN-beam. The criterion of the maximum principal stress was chosen for the initiation and propagation of the crack [8]. The direction of the crack is the direction of rotation of the stress in the element to the principal stresses. When the load grows, the stress in the model increases, and the material effort increases with them. Material effort (μ) is a quantity described by Podgórski in [9]. If the effort exceeds 1, then the crack goes through the next finite element. The calculated material effort depends on the chosen element type. For X-FEM in Abaqus, only four-node elements are available. When the reduced integration is chosen, then the material effort is calculated as the ratio of the maximum principal stresses in the centre of the element to the tensile strength. When no reduced integration is selected, the program calculates this effort for four integration points of the element and then determines the arithmetic mean. The same goes for determination of the direction of the crack. This is not an ideal approach, because the most correct would be to take stresses from the point of the crack tip instead of from the centre of the element. In this way more force is needed to destroy the element. It also depends on the size of the mesh. The larger the finite element, the further away from the crack tip are the integration points, and hence the smaller stresses in them. This means that for correct calculations a dense mesh is needed.

Results of the crack path and maximum force
In previous works, the authors described that the path calculated in Abaqus using XFEM is unrealistic [10]. However, in the case of the SEN-beam test, the crack path results are satisfactory, and the only drawback is that the crack cannot pass through the last finite element. The Fig. 3a presented below is the picture from the program and the Fig. 3b is the comparison with the results of laboratory tests obtained by de Borst [2] and Pan et al. [3]. a) b) Fig. 3. The crack path for SEN-beam. a) Results from Abaqus. Source: the authors' own study, b) the results from laboratory tests. Source: [3] However, the obtained maximum force is very different from the similar tests presented in the source literature [3], [5]. The maximum force in the simulation is 110.1 kN, and for laboratory tests, the force is about 40 kN.

The main algorithm
Abaqus allows implementation of the custom crack propagation criterion with the use of user subroutines. The authors prepared several subroutines in Fortran. These are two algorithms, the first concerns the own approach to the criterion of maximum principal stresses, the second is the implementation of the Ottosen-Podgórski criterion [9]. The main simplified algorithm of Abaqus operation with the use of subroutines used by the authors is shown in Fig. 4. First, Abaqus calculates the model and then uses the UVARM subroutine, and then the UDMGINI subroutine separately for every integration point and every iteration. The UVARM subroutine is needed only to read the stresses in every integration point and transfer them by "Fortran Common block" to the UDMGINI subroutine. This subroutine determines the crack direction and the value of the searched criterion. When it decides to lead the crack to the next element (solution found), it writes the results in the Result File. When the URDFIL subroutine reads the Result File and all calculations for determining the crack are performed in this subroutine, it transfers the obtained crack propagation angle to the next increment with the Common block [8]. The reason that all calculations for determining the crack angle are made in URDFIL subroutine, not in UDMGINI is that the URDFIL subroutine is performed only every increment which is dozens of times less than UDMGINI subroutine, which is performed for every integration point and during several iteration for one analysed finite element. Also, the URDFIL subroutine is performed once for the whole model, not for enriched finite elements. This procedure is possible because within one increment the stresses change in relation to each other, and the angle is determined by proportions of the stresses in the model, not by their values, which is explained below.
As mentioned earlier in the case of four integration points in the element, the criterion is averaged from these four points. In the described method, in the first three integration points Abaqus only reads the data in the UVARM subroutine needed for calculations, and in the UDMGINI subroutine the effort is set to 0. However, all calculations are made only at the fourth integration point, and then the most efforted integration point is selected and the effort value is set 4 times bigger; therefore, after calculating the average by Abaqus, the effort value obtains the desired value. This can be illustrated by the example of the following formula:

Xxxx -to be completed during the formatting process
When the URDFIL subroutine reads the Result File and all calculations for determining the crack are performed in this subroutine, it transfers the obtained crack propagation angle to the next increment with the Common block [8]. The reason that all calculations for determining the crack angle are made in URDFIL subroutine, not in UDMGINI is that the URDFIL subroutine is performed only every increment which is dozens of times less than UDMGINI subroutine, which is performed for every integration point and during several iteration steps for one analysed finite element. Also, the URDFIL subroutine is performed once for the whole model, not for enriched finite elements. This procedure is possible because within one increment the stresses change in relation to each other, and the angle is determined by proportions of the stresses in the model, not by their values, which is explained below.
As mentioned earlier in the case of four integration points in the element, the criterion is averaged from these four points. In the described method, in the first three integration points Abaqus only reads the data in the UVARM subroutine needed for calculations, and in the UDMGINI subroutine the effort is set to 0. However, all calculations are made only at the fourth integration point, and then the most efforted integration point is selected and the effort value is set 4 times bigger; therefore, after calculating the average by Abaqus, the effort value obtains the desired value. This can be illustrated by the example of the following formula: The same thing applies to determination of the angle of the crack. There is no simpler solution, because Abaqus does not allow to modify its built-in algorithm of separate calculation for all integration points. The method described above is not ideal because it takes stress at integration points. Ideally, stress should be applied exactly at the crack tip.

Determining the crack angle
To determine the direction of the crack, the authors used Griffith's crack and Westergaard solution [11]. Griffith's crack is shown in Fig. 5. Griffith's crack. Source: [11] Stresses around the crack tip depend on material parameters, angle θ (from -180° to 180°), distance from the crack tip r, and can be simplified to: The same thing applies to determination of the angle of the crack. There is no simpler solution, because Abaqus does not allow to modify its built-in algorithm of separate calculation for all integration points. The method described above is not ideal because it takes stress at integration points. Ideally, stress should be applied exactly at the crack tip.

Determining the crack angle
To determine the direction of the crack, the authors used Griffith's crack and Westergaard solution [11]. Griffith's crack is shown in Fig. 5. Fig. 5. Griffith's crack. Source: [11] Stresses around the crack tip depend on material parameters, angle θ (from -180° to 180°), distance from the crack tip r, and can be simplified to: After calculating the material effort for principal stress criterion and the Ottosen-Podgórski criterion [9] for r = 1, the graph is obtained (shown in Fig. 6). As it is a symmetrical task, it can be concluded that the crack will be led at an angle where the minimum material effort value is obtained. Fig. 6. Relation between the material effort for both criteria and angle relative to the crack. Source: the authors' own study based on [9].
In Abaqus, it is not possible to read the crack tip coordinates directly from the Results file. Therefore, before the presentation of the operation of the crack propagation algorithm, the data from Abaqus PHILSM should be explained. PHILSM is a result that occurs only on cracked finite elements. It has positive values on the one side of the crack and negative on the other. With a linear fit, the crack is where the PHILSM takes the value 0. Knowing the PHILSM values in the nodes and their coordinates, the algorithm determines the coordinates of the crack tip. As it was said before, the crack angle is calculated in UVARM subroutine. The outline of its process is as follows:  Reading the "Result file" for the previous increment for the whole modelcoordinates and stresses in integration points, PHILSM, and coordinates for nodes;  Calculating the material effort for all integration points. Material effort formulas depend on the chosen destruction criterion. The value of material effort depends only on stress;  Calculating the crack tip coordinates from the PHILSM values;  Reducing material effort values only to about 30 nearest integration points around the crack tip;  Bringing material effort values to the radius = 1, knowing that: 5 Fig. 5. Griffith's crack. Source: [11] Stresses around the crack tip depend on material parameters, angle θ (from -180° to 180°), distance from the crack tip r, and can be simplified to: (2) After calculating the material effort for principal stress criterion and the Ottosen-Podgórski criterion [9] for r = 1, the graph is obtained (shown in Fig. 6). As it is a symmetrical task, it can be concluded that the crack will be led at an angle where the minimum material effort value is obtained. Fig. 6. Relation between the material effort for both criteria and angle relative to the crack. Source: the authors' own study based on [9] In Abaqus, it is not possible to read the crack tip coordinates directly from the Results file. Therefore, before the presentation of the operation of the crack propagation algorithm, the data from Abaqus PHILSM should be explained. PHILSM is a result that occurs only on cracked finite elements. It has positive values on the one side of the crack and negative on the other. With a linear fit, the crack is where the PHILSM takes the value 0. Knowing the PHILSM values in the nodes and their coordinates, the algorithm determines the coordinates of the crack tip. As it was said before, the crack angle is calculated in UVARM subroutine. The outline of its process is as follows: • Reading the "Result file" for the previous increment for the whole model -coordinates and stresses in integration points, PHILSM, and coordinates for nodes; • Calculating the material effort for all integration points. Material effort formulas depend on the chosen destruction criterion. The value of material effort depends only on stress; • Calculating the crack tip coordinates from the PHILSM values; • Reducing material effort values only to about 30 nearest integration points around the crack tip; • Bringing material effort values to the radius = 1, knowing that: 6 Fig. 7. Explanation of the PHILSM values. Source: the authors' own study.
As it was said before, the crack angle is calculated in UVARM subroutine. The outline of its process is as follows:  Reading the "Result file" for the previous increment for the whole modelcoordinates and stresses in integration points, PHILSM, and coordinates for nodes;  Calculating the material effort for all integration points. Material effort formulas depend on the chosen destruction criterion. The value of material effort depends only on stress;  Calculating the crack tip coordinates from the PHILSM values;  Reducing material effort values only to about 30 nearest integration points around the crack tip;  Bringing material effort values to the radius = 1, knowing that: to get the result for radius r = 1 material effort is multiplied by Xxxx -to be completed during the formatting process to get the result for radius r = 1 material effort is multiplied by √ ,  Some values on the obtained graph are shifted 360 degrees to the left or right so that the angle from the previous increment is 0;  Fitting the 6th-degree polynomial to points using the method of least squares;  Finding the minimum value nearest to the angle from the previous increment by the bisection method. The new angle of the crack tip is the x value where this minimum was found. The last four points listed above are illustrated on exemplary values in Fig. 8.

Results for the authors' own subroutine
The above subroutine was submitted to the same simulation as in chapter 3.2. For the same criterion (maximum principal stress), the destructive force is 108.7 kN, and for the Ottosen-Podgórski criterion, the force is 104.9 kN. For both criteria, the crack paths almost do not differ from each other (Fig. 9a). The Fig. 9b shows the differences between crack path for default Abaqus criterion, for the result with using subroutines, and for laboratory tests.
The obtained results are visibly similar, but the simulation with user subroutine gives outputs more similar to laboratory tests. In addition, some correlation was also noticed -the further to the right the crack reaches, the more strength is needed to destroy the model. If the crack obtained by the default Abaqus method were closer to the laboratory results (more to the right), the value of 110.1 kN would be even greater, and thus even further away from laboratory tests. This means that the results obtained with the use of the subroutine are closer to reality both in terms of the crack path and the maximum force obtained.

Results for the authors' own subroutine
The above subroutine was submitted to the same simulation as in chapter 3.2. For the same criterion (maximum principal stress), the destructive force is 108.7 kN, and for the Ottosen-Podgórski criterion, the force is 104.9 kN. For both criteria, the crack paths almost do not differ from each other (Fig. 9a). The Fig. 9b shows the differences between crack path for default Abaqus criterion, for the result with using subroutines, and for laboratory tests.
The obtained results are visibly similar, but the simulation with user subroutine gives outputs more similar to laboratory tests. In addition, some correlation was also noticed -the further to the right the crack reaches, the more strength is needed to destroy the model. If the crack obtained by the default Abaqus method were closer to the laboratory results (more to the right), the value of 110.1 kN would be even greater, and thus even further away from laboratory tests. This means that the results obtained with the use of the subroutine are closer to reality both in terms of the crack path and the maximum force obtained. a) b) Fig. 9. a) A crack path for SEN-beam obtained with the use of the subroutine. Source: the authors' own study, b) The comparison of all crack paths: Source: the authors' own study based on [3] Despite the above conclusions, the obtained force results are not satisfactory. Although they are better than the results for the default method, they are definitely different from the laboratory results.
However, it is worth paying attention to the diagram in Fig. 10. It presents the relation between the degree of cracking and strength during calculations. The crack length is expressed as a percentage ratio of the current crack length to the final length (total failure). This figure shows the cyclical increase in force and crack as well as the decrease in strength and stop of the increase in crack length. It should be noted, however, that at the end of the second increase in force, when the force was about 55 kN (close to the force from laboratory tests), the crack reached 90% of its final length, which is almost full destruction. Then, with each large increase in strength, the crack almost did not grow. The final part that should break is the area where very low tensile forces act on the crack, and compressive forces dominate, which are not taken into account in the discussed failure criteria. Another reason for this phenomenon may be a different type of support in comparison with the performed laboratory tests. Perhaps it might be influenced by the ability of the lower supports to move vertically. In the future, the authors plan to change the boundary conditions and replace them with additional elements in the model, that allows all movements occurring in the laboratory test. Simulation of four-point beam bending test using the X-FEM method 61

Summary
The described calculations show that computer analysis gives results that are far from reality in terms of destructive force. SEN-beam is an interesting laboratory test in which many factors influence the results. It is particularly important to examine what phenomena occur in the final phase of the test when the crack is near the opposite side. For this particular test, the results obtained from tests made with user subroutine are slightly better than the ones with the use of the default maximal principal stress criterion.
It is also interesting that for this particular experiment the results obtained by default in the Abaqus program are very similar to those obtained with the subroutine. Usually, the calculations using this subroutine give much better results. For example, in a simple three-point bending test of a notched beam shown in Fig. 11, it is clearly visible that the results obtained from default in-build Abaqus criterion (a) are completely unrealistic, and the results from the simulation using the subroutine (b) are very good. Moreover, for the above example, the force obtained with the authors' subroutine is very close to real value in comparison with the default Abaqus criterion. But this is the topic for future papers.