The Computational Method of Predicting the Crack Path in Confrontation with Laboratory Tests

The paper presents the results of laboratory tests and computer simulations of the pull-out test of a rod embedded in a concrete cylindrical sample. The main purpose was to verify the force-displacement relationship and the crack paths. The X-FEM method of fracture simulation was used. The authors compared a simple maximum principal stress criterion, as well as their own method of simulating fracture propagation in brittle materials, implemented with the user subroutine in the Simulia Abaqus FEA system. A number of conclusions were obtained in the analyzes, related to the influence of the finite element mesh size, the method of modeling the rod and its contact with concrete. The advantages of own method were presented. the reasons for the differences between the results of simulations and laboratory tests were presented.


INTRODUCTION
The authors of this paper developed a method for simulating fracture propagation in brittle materials, implemented with the user subroutine in the Simulia Abaqus FEA system [14]. The developed method enables the implementation of any crack propagation criterion by an appropriate modification of the program code.
This subroutine acts as a criterion for determining the direction of crack propagation and cooperates with the X-FEM fracture simulation method. The X-FEM method is one of the most popular numerical techniques based on the generalized finite element method (G-FEM) and the partition of unity method (PUM). It is used to simulate a discontinued mechanical problem, especially to simulate crack growth [7]. Separation is simulated by adding an enrichment function to a finite element shape function, which is responsible for introducing discontinuities in the geometry of the finite element. The X-FEM method, therefore, allows the fracture line to be traced independently of the mesh, without using the remeshing technique.
This feature of the X-FEM method contributes to its obvious advantage of eliminating the modification of the mesh with each crack growth (remeshing), or creating models with a highly densified mesh, as in methods that use element deletion or separation. In addition, calculations using this method are relatively fast, thanks to the enrichment of the shape function of only those elements where a fracture may occur. Simulia Abaqus is by far the most popular program that can simulate fracture using the X-FEM method. The X-FEM method was implemented there in 2009 by Giner et al. [3]. At that time, it was created with the use of the UEL (user element) subroutine, which allows for adding a custom element, and now is permanently integrated with the program. Abaqus provides tools that allow the user to modify a part of the solver. There are 61 types of subroutines, including the above-mentioned UEL, as well as the most popular UMAT (user material), which allows to The Computational Method of Predicting the Crack Path in Confrontation with Laboratory Tests define a custom material. The author of this article used three types of subroutines: UDMGINI (user damage initiation), which allows to define of own crack propagation criterion, and two auxiliary ones: UVARM (user subroutine to generate element output), which allows saving the results from the entire model to a file and read from this file and subroutine URDFIL (user subroutine to read the result file). The operation of these subroutines will not be described in detail in this paper. More about them can be read in the author's previous publications [4,6].
The most important subroutine here is UD-MGINI, which in each load increment reads the values of stresses, displacements, and other basic parameters from the model, and using code written in Fortran, calculates and passes to the solver the crack propagation direction and the answer if the crack propagation condition has been met. There are few papers in the literature describing own criteria for the crack propagation created by this subroutine. They mainly concern cracks in steel [8,10], composites [1,11], or polymers [2]. There are no papers describing the fracture of brittle materials, and all the above examples relate to the simple implementation of the developed crack propagation criteria. The method developed by the authors does not use any new criteria for crack propagation, but significantly improves the method built into the Abaqus system.

DESCRIPTION OF THE PROBLEM
So far, the authors of the paper have checked the effectiveness of the described method in the three-point bending test of a notched beam and the pull-out test of a self-undercutting anchor [6], a four-point bending test of a notched beam [4] and the simulation of the Griffith's crack [5]. The main goal of this paper is to check how the developed method deals with the determination of the crack path in the task presented below. It is a concrete cylinder with a nominal diameter of 150 mm and a height of 300 mm, with symmetrically embedded steel threaded rods on both sides, along the height in both cylinder bases. The test consists in applying an axial tensile force at the end of both rods, which will result in the rods being pulled out with a fragment of concrete. The scheme of this task is presented in Figure 1. For the purposes of this work, laboratory tests of the presented study and computer simulations were carried out using the default fracture propagation criterion built into the Abaqus system and using own method to compare the pull-out force and the shape of the crack path.
There are 3 different anchoring lengths: 10 cm, 7.5 cm, 5 cm, and 2 types of anchoringthe threaded steel rod itself (variant called "T" for simplicity) and a threaded rod with a nut attached at its end ("N" variant). In the first case, the force from the rod will be transferred to the concrete through the adhesion of the concrete to the thread, because the rod is placed in the sample before the concrete setting stage. In the second case, the force is transmitted partly through the adhesion of the concrete to the thread, but mainly through the vertical pressure of the nut on the concrete. It should also be noted that the anchorage depth in the first case is measured from the surface of the cylinder base to the end of the rod, while in the second case to the beginning of the nut (where the nut presses against the concrete).
One of the initial problems was to select the diameter of the rod and the steel grade of the rod so that the steel would not plasticize during the laboratory test. The problem was, of course, the unfamiliarity with the tensile strength of the concrete before making all the samples, therefore, for safety reasons, concrete C45/55 was adopted. The average tensile strength is f t = 3.5 MPa, which, when multiplied by the cross-section area of the cylinder, gives the force F = 61.85 kN. The minimum rod core diameter is determined from the formula below: where: k r -the allowable stress in the rod, i.e. the yield strength of the steel.
The steel grade of the bolt was assumed to be 8.8, for which k r = 640 MPa. For the above data, d 1min = 11.109 mm, which is met by the M14 thread [15] with a small margin. For safety, the M16 rod was adopted. Additionally, the load capacity of the thread was checked for samples with a rod and a nut. The minimum height of the nut is: where: p is the thread pitch (2 mm for an M16 thread), d -bolt diameter (16 mm), D 1 -bolt core diameter (13.835 mm), k d -rod yield point (320 MPa for a typical 4.8 grade steel). For the above data, H min is 7.62 mm and the height of the M16 nut is 13 mm, which means that no plasticization of the thread will occur.

Basic mechanical parameters
To ensure that the simulation is as similar to the laboratory test as possible, in addition to the same geometry, the material parameters are also important. Therefore, tests of the concrete from which the samples were made were planned. Concrete was adopted from the following recipe: • cement -350 kg/m 3 , • sand -676 kg/m 3 , • aggregate 2-8mm -1205 kg/m 3 , • plasticizer -3 kg/m 3 , • water -150 l/m 3 .
It is planned to use an elastic, isotropic, homogeneous concrete model in the simulations. This means that the material parameters necessary to simulate cracking in the Abaqus system are compressive strength, tensile strength, Young's modulus, Poisson's ratio, and critical strain energy release rate in mode I. The tests were performed 28 days after the samples were made.
First, the compressive strength of concrete was measured during a compression test of 5 cubes with nominal dimensions of 150 × 150 × 150 mm. The average compressive strength of f c = 44.4 MPa was obtained with a 2% ratio of standard deviation to average.
The tensile strength was tested using the "Brazilian" method. Five cubes with nominal dimensions of 150 × 150 × 150 mm were subjected to the tensile splitting test. The average tensile strength of this test was obtained as for cylindrical specimens from the following formula: where: P max is the maximum force, d -cylinder diameter (in this case cube height), h -cylinder height (also cube height). The obtained average tensile strength was 3.273 MPa with a standard deviation to average ratio of less than 1%.
The Young's modulus was determined after performing pull-out tests on three sections of cylindrical samples with dimensions of 150 × 150 mm. To determine it, a modimeter equipped with two extensometers was used. The tests were carried out in accordance with the ASTM C 469-02: 2004 standard [13], and the module was calculated from the formula: where: ΔP -the increase in force in the last stage of the test, A -cylinder cross-section, Δu -increase in displacement in the last stage of the test, averaged over two extensometers, L -base length for the modimeter (80 mm). The mean value of Young's modulus of 38.206 MPa was obtained with the standard deviation to mean ratio of 0.5%.
Critical fracture energy was estimated based on international guidelines CEB -FIB Model Code 1990 [12] according to the formula: where: f cm -the average compresive strength, f cm0 = 10 MPa is base strength and G F0 is a base value of fracture energy depends on the maximum size of the aggregate grains, which for 8 mm aggregate is G F0 = 0.025 N/mm. For these data and assumption of f cm = 44.4 MPa the estimated critical fracture energy value is G F = 0.071 N/mm.
The Poisson's ratio was adopted, typical for concrete, ν = 0.2. Photographs of the tests are shown in Figure 2.

Pull-out test
The main study carried out for the purpose of this work was the pull-out test ( fig. 3a), which was then simulated using the finite element method. Three samples were made for each type of anchorage (T -only rod, N -a rod with nut) and for each anchorage length (5 cm, 7.5 cm, 10 cm), which gives 18 samples in total.
The concrete samples were made in typical cylindrical forms, but plates with threaded holes were used to accommodate the rods. As this is assumed to be an axially symmetrical task, the most precise alignment of the rods had to be ensured. Each displacement of the rod from the axis or its rotation could cause shear forces during pulling out, which could lead to uneven tearing of the concrete fragment, and the force needed for pulling out could be lower than the theoretical force with perfect alignment. Hence the need to make a thread in the base of the form.
Additionally, to further reduce the eccentric effect of the rods, it was decided to use an unusual joint on one side of the sample. The model of the joint is shown in fig. 3b. It was assumed that the loose in 10 nuts with washers would be sufficient to fulfill the role of a joint. The joint consists of 5 M16 threaded rods, nuts, washers, and two 1 cm thick steel plates made of S350 steel. According to preliminary computer simulations, it was estimated that the steel in the joint should not become plasticized during the pull-out test.
The load was controlled by a vertical displacement. Changes in force and displacement were monitored during the test. It turned out that the test is proceeding very rapidly. The force increases approximately linearly, followed by a brief flattening of the force diagram, and finally, the destruction of the sample occurs with a rapid drop in force. With a certain force, the crack initiation occurs at the beginning of the anchorage length. As the crack grows, the area remaining to break decreases (the area of the stress field holding the concrete fragment against breaking out), but the distance between the axis of force and this stress field increases. A short flattening of the force diagram suggests that as the crack length increases, the product of the stress volume and the moment of force is not the largest at the very beginning. Only when a certain crack length is reached it leads to a full rapid fracture because the maximum moment of force has been reached. Detailed graphs of the relationship between the force and vertical displacements will be presented later along with the results of computer simulations. After tests were performed, it turned out that a significant part of the pulled-out fragments was symmetrical, which indicated a fairly accurate alignment of the rods during the forming of the samples (example in Fig. 4a). However, part of the fracture surfaces was formed asymmetrically (example in Fig. 4b), which indicates the eccentricity of the rods, or is caused by inhomogeneity in the concrete.
Two interesting phenomena were also observed. In the case of the specimens with the shortest anchorage length, additional splitting of the pulled-out fragment often appeared in the form of three vertical cracks, usually into 3 equal parts (see Fig. 5a). The reason for this is obvious. As the rod is pulled out, the concrete at the point where the force is applied moves upwards, which causes the concrete on the top surface to move horizontally from the rod. This displacement causes tension around the perimeter of the concrete in the vicinity of the rod. In the case of short anchorage lengths, the vertical fracture area (vertical cross-section of the pulled-out fragment) is small, so the tensile force causes a vertical fracture. In the case of longer anchorage lengths, the vertical cross-sectional area of the pulled-out fragment is large enough to prevent this splitting from occurring. Nevertheless, this splitting most likely occurs at the moment of a sharp drop in force, so no significant differences in the results were noticed between the samples that were split and those that did not.
The second phenomenon, unfortunately, significantly complicates the analysis of the discussed task. It is shown in Fig. 5b. For some samples, branching of the fracture appeared near the side surface of the cylinder. This is especially troublesome because the X-FEM method used in the simulations does not allow the fracture to branch, therefore it is not clear which of the cracks obtained in the tests should be considered appropriate in the comparison with the simulations. This effect is visible also in Fig. 4, where it can be seen that the crack runs in one straight line and turns downward near the end.
All 18 pulled-out fragments were photographed. Based on the photos, several crack lines were drawn for each type of anchorage and anchorage length. In addition, 4 selected samples were scanned in a 3D scanner, which allowed for more accurate drawings of the fracture lines for these 4 samples based on the point cloud obtained from the scans. These scans are shown in Fig. 6.

ABAQUS USER SUBROUTINE Abaqus built-in method
Before the results of computer simulations are described, the authors decided that it would be necessary to describe the `operation of the builtin criterion in Abaqus and then own method of predicting the direction of the fracture.
The sequence of the Abaqus solver steps will now be briefly described. The simulation is divided into steps defined by users. In the case of a simple task, as described in this paper, it is only one step. The steps are divided into increments. This is the entire calculation step divided into force increment stages, which allows, for example, an accurate reading of stresses as the load increases. All simulations carried out for the purposes of this work were divided into approximately 100 increments. Each increment is divided into iterations. These are all attempts to find a solution in one increment. In the described simulations, 1 to 6 iterations per increment were required.
The default criterion that works with the X-FEM method is the maximum principal stress criterion. It works in a simple way. Abaqus reads the stresses at all integration points in the elements enriched for the X-FEM method in each iteration. These are the 3 elements closest to the crack tip (see fig. 7). The four-node element CPS4 has 4 integration points, while the reduced four-node element CPS4R has 1 integration point. For an axially symmetric task, these are the elements CPX4 and CPX4R, respectively.
Then the stresses at the integration points are rotated to the principal stresses. The angle of rotation of the stress tensor α is also the angle over which the crack will be led. For elements with four integration points, the angle is averaged from four values. Abaqus leads the fracture to the next finite element when in the last iteration of the increment the maximum principal stresses averaged from the four integration points exceed the tensile strength value. This is a very simple method and has several drawbacks. Only the closest integration points where the highest disturbances occur are taken into account for the crack direction calculations. In addition, stresses occurring in close proximity to the crack tip are not considered. For this reason, the crack paths determined by this criterion are often incorrect.

Author's method
The author's method is described in detail in a previous publications [5,6]. Only its assumptions will be presented here. The maximum principal stresses criterion was used in the author's method but considered differently than in the default method.
This method is based on the Westergaard solution of the Griffith's crack [9].
where: K I -critical stress intensity factor in mode I, r -distance from the crack tip, θ -the angle between the crack tip and a given point where the stresses are calculated.
These are the stresses around the crack tip in the infinite domain stretched by a force perpendicular to the fracture (see Fig. 8a). After determining the principal stresses from the above formulas, a graph is obtained as in Fig. 8b. These are the principal stresses around the crack tip at an equal distance from the crack tip. The local minimum for an angle of 0° can be seen here. The authors assumed that in the Abaqus simulations also the fracture direction would be guided toward the local minimum for any task.
To save computation time, the main code searching the crack propagation direction is implemented in the URDFIL subroutine, which runs only once per increment, and then the results are only transferred to the UDMGINI subroutine, which runs in each iteration. At the beginning of each iteration, the stresses in each integration point are read from the Abaqus "Results" file. Then the coordinates of the crack tip are calculated and all integration points are discarded, except for several dozen points closest to the crack tip. The maximum principal stress is then calculated at these points. Then, these stresses are interpolated as if they were at the same distance from the crack tip, on the basis of the same relationship as in the Westergaard solution, i.e. stresses decrease with increasing distance from the crack tip at the rate of 1/√ .
The stress results obtained in this way are analyzed by the dependence on the angle θ around the crack tip. On this basis, a graph as in Fig. 9 can be made. The sixth-degree polynomial is fitted to the points on the graph by the least squares method, and then by the bisection method, the local minimum closest to the angle from the previous increment is found. The angle for this local minimum is the crack propagation angle and is passed on to the Abaqus solver.
These subroutines are written in Fortran and are compiled each time the simulation is run. The described criterion has two significant advantages over the default method: stresses from more points than just those closest to the crack tip are taken into account. Moreover, the propagation angle in this method is almost unaffected by shear stresses, while in the built-in method, the propagation angle is the rotation of the stress tensor to principal stresses, which is strongly dependent on shear stresses, and those near the crack tip are subjected to very high disturbances.

Description of the numerical model
The rod pull-out tests described in chapter 3.2 were subjected to computer simulations using the built-in criterion and own method. The described task is axially symmetrical. Additionally, it is symmetrical in relation to the horizontal plane, therefore a quarter of the vertical cross-section of the cylinder was modeled with the axis of symmetry. The model includes the sample itself without the rod. The load was modeled by a vertical displacement at the point of contact between the rod and the concrete. Additionally, the horizontal displacement was blocked to simulate the adhesion of concrete to the thread.
An attempt was made to perform a simulation with a modeled rod, because in reality there is steel elasticity, so there is a non-homogenous displacement along the contact between the rod with concrete. This attempt is described in the chapter 5.4, however, none of the four proposed ways of simulating the rod gave better results than the simulation without the rod. For this reason, these simulations were discarded, and the focus was on simulations without the rod. However, the problems associated with the contact between concrete and steel are not caused by flaws in own procedure, but by the Abaqus algorithm.
There is a problem of whether there should be a horizontal displacement lock at the contact line. This is a permanent connection, however, in some tests, the concrete was separated from the rod. Such a phenomenon means that there should be no blocking of the horizontal displacement at the separation line, but also the forced vertical displacement. However, this phenomenon occurred only in some samples, and in addition it appeared only at the last moment of the tests, so it was decided to leave the horizontal displacement lock on the line of contact.
A vertical boundary condition has also been inserted on the bottom edge as this is the plane of  Fig. 10a.
The mesh of the model consists of 5 mm CPX4 four-node elements with a density of up to 1 mm in the area of the predicted crack path (Fig. 10b, Fig. 10c). The material was selected as elastic with the previously calculated Young's modulus and Poisson's ratio. The failure model was set using own subroutine or the maximum principal stresses criterion (MaxPS) with the above-determined tensile strength and fracture energy. The cracking was simulated using the X-FEM method without indicating the initial crack. Geometric nonlinearity and automatic stabilization are excluded. The calculations were made using the Abaqus/Standard solver.
After the first few simulations, it turned out that the obtained maximum force is several dozen percent greater than that obtained in laboratory tests. It was assumed that this problem may be caused by the assumed value of the critical fracture energy. Therefore, it was decided to perform all simulations for two variants of the critical fracture energy: G Ic = 0.071 N/mm, as assumed initially, and G Ic = 0.015 N/mm, i.e. the value for which the results were closest to the values obtained in the laboratory tests.

Computer simulation results -force diagrams
This chapter presents the simulation results using the default maximum principal stress criterion and own method in combination with the corresponding laboratory test results. The forcedisplacement relationship diagram and the shape of the crack path will be compared.
In the above graphs ( fig. 11, fig. 12, fig. 13), lines named "Test" are laboratory test results, "Default" are simulation results using the default method, and "Own" are simulation results obtained using the author's user subroutine. In addition, the values of displacements in laboratory tests turned out to be many times greater than those in simulations, therefore they were reduced fifty times for better readability.
As you can see, there are many differences between the simulation results and the laboratory test results. The maximum force obtained in the simulation with the critical fracture energy of 0.071 N/mm turns out to be from 33% to 75% higher than the laboratory test average (greater for a shorter anchorage length). For this reason, a lower fracture energy of 0.015 N/mm was assumed, for which the maximum value of the force differs by 1-2% from that obtained from the tests. However, it cannot be assumed that this is the actual fracture energy of the tested material, and it was only selected to match the results of laboratory tests. What's important, this value fits all test results. It is possible that other factors also influence the difference in the maximum force value. First of all, it should be mentioned that the simulation is perfectly symmetrical. In the case of laboratory tests, the rods probably did not fit on one axis. A small eccentricity could cause a shear force which can cause the rod to be pulled-out sooner than it should have been. It is also possible that the boundary conditions were incorrectly modeled in the simulations. This will be verified later in the paper by adding a rod and a contact between the rod and concrete to the simulation in various ways. Another reason for the difference may be caused by the criterion that determines the fracture. It is based on stresses that theoretically escape to infinity at the crack tip, and thus depend on the size of the mesh. Also, in the further part of the paper, the influence of the mesh size on the simulation results will be analyzed.
Another difference between the simulation results and the laboratory results is the nature of the failure. As can be seen in the charts -laboratory tests are highly dynamic. Displacement increases linearly, followed by sudden failure and a rapid decline in force. In simulations, after reaching the maximum force, it decreases slowly. As mentioned at the beginning, in reality, the force and crack path increase until a certain critical force is obtained. When it is exceeded, a quick fracture occurs, as no more force is required to guide the crack further. In simulations, although the force decreases, the displacement must continue to increase for the fracture to progress. It was suspected that this was due to too few load increments as Abaqus could not break through more than one element in one increment. Therefore, it was decided to check how the simulation would behave with a tenfold increase in the number of increments, so that in each increment a break could occur by another finite element and that, after exceeding the maximum force, the model would completely break 10 times faster. However, this did not happen. The fracture continued to progress with increasing displacement, and no significant differences were noted. Perhaps the reason is still the chosen criterion of crack propagation. When the stress exceeds the tensile strength around the crack tip, Abaqus leads it to the next element. For this reason, the crack tip moves away from the model axis and the stresses decrease a little and the program has to increase the load again for the stress at the crack tip to increase. The last difference between the results is the size of the displacements. As already mentioned, the values of displacements in the graphs in laboratory tests have been reduced fifty times, as they are about 50 times greater than the displacements in the simulations. The most likely reason for this is the way the displacements are measured during the test. The displacements in the sample were not measured -they were read from the testing machine. In fact, the displacements between the upper edge and the center of the specimen (as read from the simulation) will be much smaller than those from the machine, and this is due to several phenomena: loose in the machine, the elasticity of steel rods and elasticity of the joint, which certainly was significant. In the further part of the article, simulations with modeled rods will be presented to check the influence of steel elasticity on displacements.
There are a few more observations worth mentioning. The effect of the nut on the result is insignificant both in simulations and in laboratory tests, the maximum force is 3 to 11% greater in the cases with a nut, however, it has been noticed that the shorter the anchorage length, the greater the difference between the cases with and without the nut. The plot of the forces obtained in the simulations using the author's method is shown in Figure 14, where these slight differences can be seen in the results between the simulation with and without the nut. Furthermore, this diagram shows an almost linear relationship between the anchorage length and the maximum force.
The key question asked at the beginning of the work is "how do the results obtained using the default criterion differ from those obtained using own method". As can be seen from Figures 11 to 13 charts, these differences are small. The maximum force varies by 1-2% in each case without any particular relationship. It can therefore be concluded that this is within the error of the finite element method.

Computer simulation results -crack paths
Now the crack paths will be analyzed. In this paper, they are more important than the forces, because the author's procedure does not implement any other algorithm for determining the moment of the crack propagation than the default one. The differences in forces between the default method and the custom procedure are only due to the different shapes of the crack path. All the following crack lines in the simulations were obtained for the material with the assumed critical fracture energy G Ic = 0.071 N/mm.
An exemplary crack path obtained with the author's procedure for the 5 cm rod without the nut is shown in Figure 15. As can be seen, the default method can lead the fracture to the end, which was not often the case in the author's previous analyzes [6]. It is also visible that the fracture line is very similar to the results of the laboratory tests (Fig. 4). Here a very similar starting angle of the crack, propagating in a straight line until it is flattened near the end can be seen. Figure 16 shows a comparison of the crack paths obtained for an anchorage length of 5 cm. As you can see, the line obtained using own method has a smaller starting angle than the one obtained using the built-in criterion. Probably for this reason the maximum force for own method was usually slightly higher, as a larger fragment of the sample had to be extracted. Also, for the simulation with a nut, the crack starts a little lower for obvious reasons, but its shape itself does not depend on the presence of the nut. It was suspected that the use of the nut would completely change the nature of the load acting on the specimen, and this would affect the fracture line, but this is not the case.
The figures above (Fig. 17) show another important fact observed in computer simulations. In each variant, the fracture line flattens with increasing anchorage length, which seems logical, as if the fracture began exactly in the middle of the cylinder height, it should be completely horizontal.
At the next stage, the actual fracture lines were drawn based on the photographs of the samples after destruction. All specimens were photographed from two directions giving 4 crack paths per specimen (fracture planes every 90°), which gives 12 crack paths for one type and length of anchorage. The crack paths were drawn manually over photos, so it is not an ideal method. However, the photos were taken from a distance to reduce distortion caused by the perspective effect, and the lines themselves were drawn with proper care. As can be seen in fig. 18, the shape of the crack line itself, both for the default method and the own procedure, is very close to the actual lines. However, it can be seen that the crack path for the author's procedure is closer to the average  of the actual lines than the line obtained using the default method. So, this is a slight advantage of own method over the built-in one. Unfortunately, concrete is a heterogeneous material and the specimens themselves were relatively small, resulting in the fracture lines being very warped and having a large spread. For this reason, assessing the effectiveness of own method based on this test may be considered questionable.
In the next step, the crack lines from the simulation were compared to the lines obtained from 3D scans of individual samples to eliminate the influence of errors resulting from manual drawing and inaccuracies of photos. It should be noted, however, that the following paths were obtained from scans of individually selected samples, and these may differ from the average of three samples for one type and length of anchorage. The lines were obtained from a point cloud. On their basis cross-sections in vertical planes, every 30° were made using Blender software. This gives 12 crack paths per sample. Thus, it is a sufficient number of cross-sections to visualize the spatial surface of the pulled-out fragment based on a flat drawing.
A comparison of the crack paths is shown in Figures 19 to 22. As can be seen, the crack paths in samples 10 and 14 confirms what has been proven earlier -the line obtained with the use of own criterion is closer to the average shape of the actual fracture than the line obtained by the builtin method. Unfortunately, for sample 3, the opposite is true, but all the lines are close enough to each other that can be put within the error limits.
When it comes to sample 18, the crack path from the built-in criterion is also closer to reality, however, the initial angle of the crack in this sample is unusually large. This causes the ends of the crack line to be closer to the top plane of the cylinder, and therefore closer to the crack path for the default method.
It is also worth noting the relationship between the length of the anchorage and the height of the fracture line in the results of laboratory tests. For samples with a large anchorage length, the fracture surface is flatter (samples 3 and 10), and for shorter lengths, it is higher (samples 14 and 18).
It was also decided to compare the fracture lines obtained in the simulations with the assumed   Figure 23. Despite the significant impact of the fracture energy on the force values, it turns out that it has almost no effect on the crack paths. However, it should be noted here that Abaqus was unable to fully carry out the fracture line when it used the default crack propagation criterion. In this case, the simulations ended with an error, like most of the simulations presented in the author's previous publications [6]. So, this is another advantage of own procedure, which showed no such problem.

Other considerations about simulation
For the purposes of this article, it was decided to verify a few more conclusions and perform additional simulations. In this chapter, only simulations with a nut and an anchorage length of 5 cm were analyzed, moreover, the simulations were performed with the fracture energy G Ic = 0.071 N/mm.
In the beginning, it was decided to add the presence of a rod to the simulation. Before the load has been simulated by a displacement applied directly to the contact point of the rod with the concrete. From now, the load will be applied to the upper end of the rod, while different methods of connecting the rod to the concrete will be analyzed. 4 additional methods of rod modeling and connection to concrete were considered: 1. Element with Young's modulus 143.5 GPa, contact realized by friction with friction coefficient μ = 0.1, 2. The same element, but with a modeled thread shape, A rigid element with a modeled thread shape, 4. The same element, but the contact is made with a permanent connection between concrete and steel.
The first variant with this unusual Young's modulus arose from the necessity to create an element that replaces the stiffness of a joint present in a laboratory test. The value of this modulus was estimated by comparing the stresses in the actual modeled joint with Young's modulus 210 GPa and its replacement bar as in Figure 24. The necessity to use a bar with equivalent stiffness is due to the fact that in a planar axially symmetric problem it is not possible to add three-dimensional elements.
The second variant was modeled based on the actual catalog dimensions of the M16 thread [15]. In assumption, this thread would be able to correspond to the actual load transfer from the rod to the concrete through the vertical pressure of the thread. Unfortunately, it turned out that in the simulation, the tension of the steel along the length of the thread also tensed the concrete in this region, causing a crack as in Figure 25.
For the above reason, it was decided to abandon this variant and a third variant was proposed, where the rod is infinitely stiff to prevent the concrete from stretching along the length of the thread. Unfortunately, the contact simulated by friction between concrete and steel differs slightly from the contact in reality.
The fourth option is the same, with the difference that the contact is made as a permanent connection between concrete and steel to simulate the actual bond between concrete and steel. Fig. 26 shows the force-displacement relationship.
As you can see, the variant with a rigid rod does not differ much from the variant without a rod, while in the case of an elastic rod, the displacement at the end of the bar increases, but the maximum force remains around the same.
Unfortunately, it is impossible to determine the correct displacements in variant 1, because either exactly the same joint need to be modeled in the simulation or displacements using strain gauges on the sample need to be measured. In this simulation, however, a different conclusion is important -the nature of the force-displacement dependency curve, in this case, is almost identical to that in laboratory tests (except for the difference in the values of forces and displacements themselves). As in the tests, the line is constantly at the same slope, and at the very end, there is a short flattening and a rapid decrease in force. The other lines do not fit into this character. Of course, in these simulations also no significant differences were noticed between the variants using the default criterion and own procedure.
As for the crack paths, they are shown in Figure 27. As can be seen, there are no major differences between the variants, but there is still a visible relationship -a smaller slope of crack paths obtained using the author's procedure. Only in variant 1 of the default method, an elastic rod deviates downwards.
The last performed verification was the effect of the mesh size on the results. Until now, all simulations have been performed with 1 mm elements in the area of the expected fracture line. Now, simulations with 0.5 mm, 2 mm, 5 mm, and 10 mm mesh (Fig. 28) will be checked.
As can be seen in Figure 29, the size of the finite elements near the fracture line does not have much of an effect on the results. The maximum force increases with the size of the finite elements. In the case of the default method, the force for the 10 mm mesh is greater by 23% compared to the force for the 1mm mesh, and for own criterion, it is 15%. Moreover, the difference in the results for the 1 mm and 2 mm mesh is approximately  1% for both methods. These numbers suggest that the program correctly determines the moment of crack propagation, and the differences in the results between the simulations and the laboratory tests are due to other reasons mentioned earlier, i.e. probably a wrongly adopted fracture energy. As for the crack paths, there are also no major differences here (see Fig. 30). All the differences in the crack paths are due to the obvious decrease in the accuracy of the model because the crack tip with increasing load can only jump from the edge of the element to another edge, it cannot stop inside the element, this is a limitation of the X-FEM method. A quite significant advantage of the own method over the built-in criterion is the fact that in the default method it was not possible to draw the crack path to the end for large sizes of finite elements, while the own procedure did not manage only the simulation with 0.5 mm elements. The reason for this  was the number of finite elements, making the "Results" file very large. Own procedure uses this file constantly during calculations, which made the computer overload. However, this does not mean that the algorithm of own method is ineffective because the calculations using own criterion were about 3/4 times longer than the calculations using the default method, while for the model with the size of 1 mm elements it took 3 minutes and for 10 mm elements several seconds for mid-range personal computer.

CONCLUSIONS
The article presents the results of laboratory tests and simulations of the pull-out test of a rod embedded in a concrete cylinder. During the simulation, the default maximum principal stress criterion and own method were used to determine the direction and propagation of the fracture.
In previous author's works [5,6] it was proved that simulations using the own procedure give a crack line much more correct than the built-in criterion in the Abaqus system, especially in simulations where Mode I of the crack loading is dominant. In this work, an attempt was made to check the operation of the own procedure in the pull-out test of a rod embedded in concrete.
For the purpose of the work, rod pull-out tests were carried out in laboratory on 18 samples (3 anchorage depths and 2 types of rod end). In order to verify the mechanical parameters of the selected material, 5 compression tests on cubes, 5 Brazilian tests on cubes and 3 compression tests to find the Young's modulus were carried out. Other material parameters were adopted or estimated.
Authors performed 12 computer simulations (3 anchorage depths, 2 types of rod end, 2 different values of critical fracture energy), 4 additional simulations with different methods of modeling the rod, and 4 additional simulations with different finite element sizes. All simulations were performed second time using own procedure of crack line prediction.
It turned out that the determination of the crack paths is very difficult due to the ambiguous results of the laboratory tests. The small size of the samples and the heterogeneity of the concrete resulted in the fracture surfaces being very irregular, however, the lines obtained by both the default method and the author's method fit within the range obtained in the tests. The shape of the line follows the trend visible in the tests -the initial angle of the crack is very similar, and near the end, the crack flattens. Usually, however, own criterion determined the lines closer to the mean from the laboratory tests. The maximum force in the simulations is significantly higher than that obtained in the tests. The most likely reason for this is incorrectly assumed fracture energy. For a lower value of the fracture energy in the simulations, the value of the force needed to initiate the fracture decreased. Boundary conditions also turned out to be a problem. Although different methods of modeling the load transfer to concrete did not significantly affect the results, only for the simulation with an elastic rod the shape of the force-displacement relationship curve turned out to be similar to the curve obtained in the tests. Unfortunately, the values of displacements in the simulations were read differently than in the tests, therefore it was decided not to compare the displacements.
The anchorage length of the rod has an almost linear effect on the maximum force needed to pull out. In addition, the longer the anchorage, the flatter the crack path. This can be seen in both laboratory tests and simulations. However, there is no significant difference between the results with and without the nut at the end of the rod. Moreover, it was proved that the size of the finite elements does not have a significant influence on the results.
The paper aimed to evaluate the effectiveness of own method of predicting the crack path against the default criterion in the Abaqus system. It was shown that for a large part of simulations using the built-in method, the calculations did not reach the end and ended with an error, while in the case of the author's method, this problem occurred occasionally and was rather caused by the lack of full freedom in creating Abaqus user subroutines. The disadvantage of the own method is the slightly longer computation time.
Several obstacles were encountered in the work. For this reason, similar research and analyzes are planned for the future. It is planned to carry out tests using larger samples to reduce the effect of irregular fracture surfaces caused by the inhomogeneity of the concrete. More tests are planned but limited to tests without a nut and with one anchorage depth. It is also planned to perform these tests with the use of strain gauges so that the results of displacements can be compared. Additionally, it is planned to perform a fracture energy test instead of an estimation.