Constraint evaluation and effects on selected fracture parameters for single-edge notched beam under four-point bending


 This article presents the results of analytical and numerical research focused on the numerical determination of selected fracture mechanics parameters for beams containing a crack in the state of four-point bending with the dominance of a plane strain state. Based on the numerical results, the influence of the specimen geometry and material characteristic on selected fracture parameters is discussed. By analogy to the already known solutions, new hybrid formulas were presented, which allow to estimate the J integral, crack tip opening displacement, and load line displacement. In addition, the study verified the Shih formula connecting the crack tip opening displacement and J integral, taking into account the influence of in-plane constraints on the value of the proportionality coefficient denoted as d
 
 n
 in the analysis. This article also presents the development of Landes and Begley’s idea, which allows to experimentally determine the J integral as a measure of the energy causing the crack growth. The innovative element is taking into account the influence of in-plane constraints on the value of the η coefficient, which is the proportionality coefficient between the J integral and energy A. The last sections of this article are the assessment of the stress distribution in front of the crack tip and the selected measures of in- and out-of-plane constraints, which can be successfully used in the estimation of the real fracture toughness with the use of appropriate fracture criteria.


Introduction
Beams are one of the main structural elements commonly used in engineering, especially in construction practice. They are used in civil and mechanical engineering as structural elements carrying various loads. The beams themselves have also found recognition in the use of their geometry in mechanical tests in the field of resistance strain gauge, Charpy tests, or material tests in the field of fracture toughness, in which three-point bent beams are usually used, with a deliberately introduced crack (slot) [1,2]. These specimens, commonly designated as SEN(B) (single-edge notched in bending), are the basic geometry used to assess the fracture toughness by various methods [1][2][3], as well as the basic structural elements used in the idealization of complex objects, based on EPRI [4], API [5], SINTAP [6], and FITNET [7] procedures. This geometry has been discussed many times in the professional literature in the field of fracture mechanics, in terms of the assessment of geometric constraints and other quantities characteristic of problems in the field of fracture mechanics [8][9][10][11].
In many scientific papers, there are discussions about a beam with a crack bent with a pure moment [12][13][14]. In ref. [15], the limit loads for a beam with a crack bent at the place of the crack with a pure moment were assessed, by conducting computer simulations based on the geometry of a four-point bent beam, marked in the professional literature as SEN(4PB) - Figure 1.
Simulation studies presented in ref. [15] led to the development of formulas that allow to estimate the limit load without the need to perform numerical calculations depending on the level of the yield point and the crack length. The results presented in ref. [15] were obtained by analyzing the curves presenting the force diagram as a function of the load line displacement and by assessing the plasticity of the noncracked ligament of the specimen. The SEN(4PB) geometry can be used for modeling and the analysis of fracture processes in structural elements in which pure moment bending is observed. As shown in Figure 1, between the load rollers whose distance is 2 W (where W is the width of the specimen), the shear force value is zero and the bending moment has a constant value of M g = P·W. However, the professional literature does not provide detailed references to the analysis of the fracture mechanics parameters or the aforementioned geometric constraints, which should be understood as the limitations that the material placed on the plastic strains develop under the influence of the external load. A comprehensive analysis of the SEN(4PB) beam in the field of fracture mechanics should include: • Assessment of changes in the J integral, crack tip opening displacement, and load line displacement as a function of the external load, normalized by the limit load (similar to other basic geometries in EPRI procedures [4]), with a proposal of hybrid solutions for the SEN(B) specimens given in ref. [10]. • Mutual relation between the J integral and the crack tip opening displacement using the Shih formula [16]: where δ T is the crack tip opening displacement; J is the J integral; σ 0 is the yield point; d n is the proportionality coefficient, depending on the parameters of the Ramberg-Osgood curve, yield point, Young's modulus, and stress distribution near the crack tip determined according to the HRR solutionthe d n parameter values can be determined using the computer program presented in ref. [17].
• Verification of the parameter η, which is the proportionality coefficient between the J integral and the energy necessary to estimate it based on the Landes and Begley's formula [18]: where b is the noncracked section of the specimen (b = W − a), B is the specimen thickness, A is the energy calculated as the area under the force P curve as a function of the load line displacement v LL (or the crack mouth opening displacement v CMOD ).
• Analysis of stress fields in front of the crack tip.
• Assessment of the level of selected measures of geometric constraints: ○ Q stresses being the difference between the real distribution estimated with the use of the finite element method (FEM) and the HRR solution [19,20]; ○ Mean stresses σ m , effective stresses σ eff , and stress triaxiality coefficients in the form of the parameter T z [21][22][23] or the quotient of effective stresses and mean stresses.
The assessment of the aforementioned quantities characteristic for various problems in the field of fracture mechanics requires complex numerical FEM calculations (which can be carried out for the dominance of a plane stress state or a plane strain state, or for three-dimensional problems), together with a complex analysis of the obtained results. Due to the fact that the fracture toughness, in accordance with the standard [1], is determined for the dominance of a plane strain state, for which hybrid solutions are given in refs [4,10], this article deals only with the domination of a plane strain state, for which, according to the O'Dowd and Shih theory, the Q stresses [24,25] are also determined.  Figure 1: (a) Model of a four-point bending beam with a crack with a length of a -SEN(4PB) specimen (Wspecimen width, Bspecimen thickness, 2 Ltotal specimen length, usually amounting to 4.25·W, Pspecimen loading force) [15]. (b) Distribution of load, shear forces, and bending moment for a four-point bending beam.
carried out with the use of finite element method (FEM), carrying out numerical calculations in the ADINA SYSTEM package [26,27]. The analysis was carried out for plane strain state dominance based on the developed parameterized models of bending beams. When building the numerical model, the guidelines provided by the authors of refs [28,29] were used. Due to the symmetry, only half of the beam with the crack was modeled, with appropriate boundary conditions applied in the right place.
To reflect the actual loading of the beams in laboratory conditions as much as possible, it was decided to implement the load on the beams and their support with the use of rollers (pins, supports, etc.), which are required for solving the contact issue. This means that in the numerical model, appropriate contact surfaces and appropriate groups of finite elements have been defined so that both the support of the beams and their load reflect the actual behavior of the material during experimental tests. In the case of the considered SEN(4PB) specimens, the loading roller was modeled as a half-arc, with a diameter of ϕ16 mm and divided into 90 equal twonode contact finite elements (FEs). Displacement, which was linearly increasing in time, was applied to the defined contact surface carrying out the load. The support of the SEN(4PB) specimen in the form of a pin (support) was modeled as a half of the diameter arch equal to ϕ16 mm, which was divided into 90 equal two-node contact FEs. (This gave 91 nodes on the contact surface, similar to the load-carrying roller.) The crack tip was modeled as a quarter of an arc with a radius r w within (1 ÷ 5) μm. This means that the radius of the crack tip was 40,000 and 8,000 smaller than the specimen width in extreme cases. This crack tip was divided into 12 parts with the compaction of the finite elements toward the edges of the surface (the boundary elements, depending on the model, were (5 ÷ 20) times smaller than the largest elements located in the central part of the arch). The size of the fillet radius of the crack tip was determined by the level of the external load and by the crack length of the analyzed case of the specimen. For each specimen, the apical area with a radius of approximately (1.0 ÷ 5.0) mm was divided into (36 ÷ 50) finite elements, the smallest of which at the crack tip was (20 ÷ 50) times smaller than the last one. This meant that in extreme cases, the smallest finite element located just at the top of the crack was approximately 1/3024 or 1/10202 of the specimen width W, and the largest modeling apical area was approximately 1/151 or 1/240 of the specimen width. The parameters of the numerical model were strictly dependent on the analyzed geometry (relative crack length), material characteristics, and external load.
The analysis was carried out with the assumption of small deformations and small displacements [28][29][30], and the finite element model was filled with nine-node "2-D SOLID plane strain" finite elements with "mixed" interpolation with nine numerical integration points [26,27]. An example of a numerical model used for calculation is shown in Figure 2.
In the course of numerical calculations, a constant width of the specimens was assumed W = 40 mm, and all other external dimensions of the specimens were related to their width. The specimens with four relative crack lengths a/W = {0.05, 0.20, 0.50, 0.70} were analyzed. In accordance with the recommendations of the authors of the ADINA SYSTEM package [26,27], when carrying out the analysis for the dominance of the plane strain state, the thickness B = 1 m was assumed. In the FEM analysis, a homogeneous, isotropic model of an elastic-plastic material was used, with the Huber-Misses-Hencky plasticity condition, described by the following relationship: where σ is the stress, ε is the strain, σ 0 is the yield point, and ε 0 is the strain corresponding to the yield point, calculated as ε 0 = σ 0 /E, where E is the Young's modulus, α is the hardening constant, and n is the strain hardening exponent in the Ramberg-Osgood law.
The calculations were carried out assuming the value of the strengthening constant α = 1, Young's modulus E = 206 GPa, Poisson's coefficient ν = 0.3, four values of the yield stress σ 0 = {315, 500, 1,000, 1,500} MPa, and four values of the stain hardening exponent in R-O law n = {3.36, 5, 10, 20}. This resulted in a combination of 16 hypothetical tensile curves that can be assigned according to the mechanical properties for both ferritic steels, general-purpose structural steels, and strong and weakly hardening materials [31,32]. The full numerical analysis included 64 numerical models differing in the yield strength σ 0 , the strain hardening exponent n, and the relative crack length a/W.
To present selected results of numerical calculations, the reference point for each analyzed specimen will be the value of the limit load, which was numerically estimated for SEN(4PB) specimens in ref. [15], where concise mathematical formulas were given that allow to estimate the limit load depending on the value of the yield stress, the relative crack length, and the predominance of a plane stress or strain state. In this article, considerations are carried out for the dominance of a plane strain state, for which the limit load can be estimated using one of the two formulas. One of the formulas traditionally refers in its form to the formulas that can be found in EPRI procedures [4]. According to a simple formula, the limit load can be calculated as follows: where P 0 is the ultimate load capacity calculated in kN, B is the beam thickness given in mm, σ 0 is the yield stress given in MPa, and the function f(a/W) is given by the following form: where the matching coefficients A 1 -A 4 are, respectively, A 1 = 0.01594, A 2 = −0.01121, A 3 = −0.01577, and A 4 = 0.01228 for cases with a dominance of the plane stress and A 1 = 0.02104, A 2 = −0.01385, A 3 = −0.01971, and A 4 = 0.01476 for cases with dominance of the plane strain. The matching coefficient R 2 is equal to R 2 = 0.952 for the plane stress state and R 2 = 0.995 for the plane strain state [15]. The second formula proposed in ref. [15], approximating the set of numerically estimated limit loads, has the following form: where P 0 is the limit load given in N, B and b are the thickness and length of the uncracked ligament of the specimen, respectively (b = W − a = W − a/W·W, where W is the width of the specimen and a is the crack length) given in m, σ 0 is the yield point in MPa, and the function f (σ 0 /E, a/W) depends on the quotient of the yield stress σ 0 and Young's modulus E and the relative crack length a/W. The function f(σ 0 /E, a/W) can be determined using the Table Curve 3D program [15], approximating the curvilinear surface f(σ 0 /E, a/W) with the following equation: Figure 2: Numerical model of the SEN(4PB) specimen: (a) simplified technical drawing of the beam with a hatched fragment, which was modeled; (b) full numerical model; (c) the apical area; and (d) the shape of the crack tip [15].
where the approximation coefficients A 1 -A 10 are given independently for plane stress and plane strain states in Table 1 [15]. The use of formulas (6) and (7) to estimate the ultimate load capacity causes, that the solution to be burdened with a maximum error of 7 and 2%, for plane stress and plane strain state respectively. The average error of matching is 2.47 and 1.01% for plane stress and plane strain state respectively. For cases of geometrical and material characteristics not covered by the research program, it is recommended to use the values for the two combinations, closest to the desired one, to estimate the ultimate load capacity of the SEN(4PB) specimen [15]. For the purposes of the analysis discussed in this article, it was decided to use the limit load calculated according to formulas (4) and (5) for the dominance of the plane strain state.

Characteristics of selected fracture mechanics parameters
The main quantities assessed during the numerical calculations were the J integral (which is treated as the crack driving force), the crack tip opening displacement δ T , and the load line displacement, denoted by v LL . These values were assessed as a function of the external load P normalized by the limit load P 0 . The J integral was determined using the virtual shift method [26,27], which uses the concept of virtual crack growth to calculate the virtual energy change [26,27]. In the analysis, eight integration contours were drawn through the area encompassing all FEs in a radius of length {10, 15,20,25,30,35,40, 45} FEs around the crack tip. The contour of integration was carried out in accordance with the recommendations [28][29][30].
It should be noted that the values of the J integral obtained from the mentioned integration contours were convergent. On the other hand, the crack tip opening displacement δ T was determined after carrying out the elasto-plastic FEM calculations, using the concept proposed by Shih [16], as shown in Figure 3. The analysis of all three parameters, the J integral, the crack tip opening displacement, and the load line displacement, was carried out in terms of the dependence of these parameters on the material characteristics (expressed by the strain hardening exponent n and yield stress σ 0 ) and the geometry of the SEN(4PB) specimen, which was expressed in terms of the relative crack length a/W (Figures 4-6).
The natural conclusion is that each of these three parameters increases with the increase in the external loadin the graphsnormalized by the limit load P 0 , and the rate of changes is determined by the material characteristics and the relative length of the crack. The shorter the crack length, the greater the values of the J integral, the crack tip opening displacement δ T , and the load line displacement v LL are observed at the same level of external load (Figures 4-6a)this sentence is true for the range of changes in the relative crack length from a/ W = 0.20 to a/W = 0.70. As shown, the curve of changes of the J integral and the crack tip opening displacement as a function of increasing external load, regardless of the material characteristics, for very short cracks (a/W = 0.05), does       not fit into the observed system. It can also be seen that the lower the material hardening degree (higher the value of the strain hardening exponent n in the R-O law), the higher the values of the elastic-plastic fracture mechanics parameters are obtained at the same level of external load (Figures 4-6b). The increase in the material strength (expressed as yield point σ 0 ) is also accompanied by an increase in the value of the J integral, crack tip opening displacement δ T , and load line displacement v LL at the same level of external load (Figures 4-6c).

Hybrid solutions for estimating the basic parameters of fracture mechanics
The catalog of numerical solutions, including 192 curves, presenting the changes of the J integral, the crack tip opening displacement δ T , and the load line displacement v LL , can be used to solve engineering problems, if the real object can be identified with the model SEN(4PB) specimen. Its use in the graphic form can be a problem. In 1981, the authors of the EPRI procedures [4] proposed a hybrid approach to estimate the value of the J integral, the crack tip opening displacement δ T , and the load line displacement v LL , without the need for numerical calculations, decomposing these values into elastic and plastic parts [4], as mentioned earlier [10], which provides a simplified approach to the EPRI solution [4], resigning from the decomposition of the mentioned quantities into elastic and plastic parts. Based on the analysis presented in ref. [10], for the considered geometry, to simplify the search for the value of the J integral, the crack tip opening displacement and the load line displacement as a function of the increasing external load, material characteristics, and specimen geometry, the following empirical expressions can be proposed to calculate the aforementioned fracture parameters: small dependence of the h * function, required to estimate the selected fracture mechanics parameters based on the hybrid formulas given in formulas (8)(9)(10). The yield point has the least influence on the change of the value of the h * function for the selected material and geometric characteristics. As can be observed in the tested range of external loads, the values of the h * function strongly depend on the external loadwith increasing external load, a decrease in the value of h * is observed, with a tendency to saturate the curves h * = f(P/P 0 ), which was observed in the case of SEN(B) specimens [10], used to determine the fracture toughness in a plane strain state domination [1][2][3]. In ref. [10], the saturation values for the h * function were given depending on the material characteristics and the relative crack length, as well as the obtained test results were described with simple analytical formulas.
For the h * = f(P/P 0 ) curves presented in this article, the saturation of these curves is practically not observed, so it is not possible, as shown in ref. [10], to specify specific values of the function h * , similar to the authors of the EPRI procedures [4]. The use of the developed catalog of numerical solutions is possible with a detailed description of changes in the h * function depending on the external load. By using the scheme shown in ref. [10], changes in the functions h h , 1 ⁎ 2 ⁎ and h 3 ⁎ depending on the normalized external load can be described by the following general equations:   Analysis of the selected fracture parameters for SEN(4PB) specimen  77 where the coefficients A 1 , A 2 , A 3 , B 1 , B 2 , B 3 , C 1 , C 2 , and C 3 are determined using the Table Curve 2D program, for each curve estimated on the basis of numerical results h * = f(P/P 0 ); ( Tables 2-5). The level of matching the curves to the numerical results was determined by the determination coefficient R 2 , which determines what part of the numerically obtained data corresponds to the empirical equation proposed for their description. The average fit in the set of h h , 1 ⁎ 2 ⁎ and h 3 ⁎ curves is 0.999, 0.999, and 1.000, respectively. This may prove the correct approach to the description of numerical results with simple empirical formulas. The approximation coefficients depend on the geometry of the SEN(4PB) specimen expressed by the relative fracture length a/W, the material characteristics expressed by the yield point σ 0 , and the material hardening level expressed by the power exponent n. Attempts made to approximate changes in A 1 , A 2 , A 3 , B 1 , B 2 , B 3 , C 1 , C 2 , and C 3 coefficients, depending on the material characteristics or the SEN(4PB) specimen geometry, did not lead to a satisfactory form of simplification, which does not mean that such operations cannot be performed. It is advisable to consider the changes of these coefficients as a function of two selected parameters, and then, it is recommended to make the obtained empirical functions dependent on the third of them.
3.2 Shih relationship analysismutual relation of the crack tip opening displacement δ T and the J integral As mentioned earlier, the J integral and the crack tip opening displacement are the two basic quantities used in fracture mechanics. Both values can be used in the formulation of the fracture criteria [33] or for the assessment of the strength of structural elements containing different defects by carrying out the analysis using CDF failure diagrams (crack driving force diagrams). These values, as stated earlier, were related by Shih [16] using the formula (1). In this formula, there is a d n coefficient, which is a proportionality coefficient depending on the parameters of the Ramberg-Osgood curve, yield point, Young's modulus, and stress distribution near the crack tip determined according to the HRR solutionthe d n parameter values can be determined using the computer Table 2: Coefficients of matching equations (11)(12)(13) to the obtained numerical results for materials described with the yield point σ 0 = 315 MPa  Table 3: Coefficients of matching equations (11)(12)(13) to the obtained numerical results for materials described with the yield point σ 0 = 500 MPa  Table 4: Coefficients of matching equations (11)(12)(13) to the obtained numerical results for materials described with the yield point σ 0 = 1,000 MPa program presented in ref. [17]. However, the Shih formula does not take into account the geometry of the structural element. In ref. [34], it was shown that the value of the d n coefficient changes with the change of the relative fracture length, which largely determines the level of in-plane constraints, i.e., the constraints posed by the material of the structural element, due to the development of plastic deformations under the influence of external load. Figure 10 shows example diagrams of the crack tip opening displacement denoted as δ T as a function of the J integral.
The evaluation of the numerical results presented in Figure 10 leads to the well-known conclusions. The crack tip opening displacement is directly proportional to the J integralthe proportionality coefficient is the d n coefficient. The shorter the fracture, the greater the crack tip opening displacement value. The more the material strengthens, the smaller the crack tip opening displacement value, which also decreases with the increasing yield strength. The linear relationship between the crack tip opening displacement and the J integral allows the use of a simple method to estimate the value of the d n coefficient for all considered cases. By dividing equation (1) on both sides by the length of the noncracked ligament of the specimen denoted by b (where b = W− a), we normalize both sides of equation (1): Denoting by The proportionality coefficient d n linking the crack tip opening displacement δ T and the J integral can be calculated as follows: As can be seen, the value of the proportionality coefficient d n is equal to the tangent of the angle of the slope Figure 11 shows the influence of specimen geometry and material characteristics on the value of the d n coefficient. The presented graphs of changes ( ) = δ f J T are for illustrative purposes only. A comprehensive analysis carried out by the author of this article shows that the value of the d n coefficient strongly depends on the strain hardening level of the materialthe less hardening of the material, the lower value of the d n coefficient is Table 5: Coefficients of matching equations (11)(12)(13) to the obtained numerical results for materials described with the yield point σ 0 = 1,500 MPa observedand it does not matter whether the specimen has a short or long crack or what is the yield point of the material from which it is made. The influence of the yield point is also importantthe higher it is, the greater the value of the d n coefficient. In the case of the impact of the crack length on the value of the d n coefficient, it can be noticed that it increases with the increase of the yield point and the decrease in material hardening ( Table 6).
The analysis of Table 6 shows that the values of the d n coefficient, calculated only on the basis of the HRR field, are not sensitive to the change of the specimen geometry. Significant differences between the results of numerical calculations and the values obtained based on the parameters of the HRR field, expressed in Table 6 by the error denoted as λ, are observed for each material model in the case of very short cracks (a/W = 0.05) at an average level of about 15%. The smallest differences in the entire set of results are observed for specimens with short cracks (a/W = 0.20).
The decrease in the material hardening degree causes that the difference between the numerically determined d n coefficient values and values based on the HRR field increases, and for selected cases, it is greater than 25%. As mentioned earlier, the application of the Shih formula [16] to determine the crack tip opening displacement should use the verified values of the proportionality coefficient d n , preferably verified by numerical calculations or also by experiment. It should be borne in mind that in the case of other specimens with a predominance of bending (e.g., SEN(B) and C(T)the SEN(B) sample was analyzed in ref. [34]) or for specimens with a predominance of tension, the values of the d n coefficient will be completely different.
The numerically calculated values of the d n coefficient are visualized in Figure 12. The presented surfaces d n = f(σ 0 , a/W) for successive materials with different degrees of hardening were approximated by the following simple mathematical formula:   Analysis of the selected fracture parameters for SEN(4PB) specimen  81 where the coefficients D 1 … D 10 are presented in Table 7, together with the coefficient of determination R 2 . The numerical results with the approximating surface are shown in Figure 12. The use of formula (18) requires the knowledge of the strain hardening exponent n, the yield stress σ 0 inserted in MPa, and the relative crack length a/W. The maximum difference between the numerically determined value and that obtained on the basis of the proposed formula (18) is 0.026 for selected points for materials characterized by weak hardening (n = 20). The conducted analysis proves that geometric constraints, the in-plane constrains, not only determine the level of the J integral or the crack tip opening displacement. It should be remembered that both of these parameters can be used in the construction of fracture criteria or in assessing the strength of structural elements containing defects. The geometrical constraints depending on the geometry of the structural element and the material characteristics affect the level of stresses in front of the crack tip [8,9], as well as the values directly related to the parameters of the mechanical fields around the crack tip, as mentioned earlier.

Evaluation of the η coefficient in the Landes and Begley formula for the calculation of the J integral
The method of determining the J integral, presented almost 50 years ago by Landes and Begeley [18], presented in this article by formula (2), has received many modifications.
The method formulated by the authors [18], using many specimens in subsequent normative documents, was adapted to the scheme based on the one-specimen method, using the potential fall-off technique or the compliance change technique [1][2][3], where in refs [1,3], the analysis of the J integral is carried out taking into account the decomposition into its elastic and plastic part, while in ref. [2], this decomposition is not taken into account, and to calculate the J integral using formula (2), it is necessary to know its geometry, the energy A defined as the area under the curve of the force P plotted as a function of the load line displacement v LL , as well as the coefficient η, which is usually equal to 2 for SEN(B) specimens. Standards [1][2][3] do not provide the value of the η coefficient for SEN(4PB) specimens. These standards contain only a formula that allows to estimate its value for the compact specimensthe C(T) typeas discussed in ref. [35], for which the η coefficient is calculated as follows: where b 0 is the initial length of the noncracked section of the specimen.  The analyzes carried out by many researchers [36][37][38][39] proved that the value of this coefficient depends on the geometry and material characteristics. The authors of refs [36][37][38][39], looking for the correlation between the η coefficient and the specimen geometry and material characteristics, always decomposed the J integral into an elastic and plastic part, as shown in refs [1,3]. References [10] presents alternative equations to EPRI formulas [4] that allow to estimate the selected parameters of elastic-plastic fracture mechanics without the necessity to carry out numerical calculations, without taking into account the decomposition into elastic and plastic components given in refs [1,3]. This approach seems to be justified because many researchers decide to use the aforementioned Landes and Begeley formula in the case of brittle fracture analysis, at the expense of the procedure of determining the stress intensity factor, analyzing the area under the P = f(v LL ) curve (force versus load line displacement) [13]. Sample diagrams of changes in the force P as a function of the load line displacement v LL are shown in Figure 13.
The mutual arrangement of the P = f(v LL ) curves will not be commented on in this articledrawing obvious conclusions is left to the readerthe obtained arrangement of the P = f(v LL )curves is as expected. It should be noted that only in terms of absolute values, the highest values of the force at the same level of load line displacement are observed for the SEN(4PB) specimens with the shortest cracks ( Figure 13b) and for specimens characterized by the material with the highest yield point (Figure 14c). However, if in the graphs shown in Figure 13, on the ordinate axis, we replace the force (given in absolute values) with a relative valuethe value of the force P is normalized a) ) c ) b  by the limit load P 0 (which, as we know, is different for each specimensit depends linearly from the yield point and strongly on the crack length), the nature of the P/P 0 = f (v LL ) curves will be completely different, as shown in Figure 14. A similar analysishowever, for typical and compliant with standards [1-3] SEN(B) specimens, was presented in ref. [13].
The inverted nature of the P/P 0 = f(v LL ) curves is visible during the analysis of Figure 14a, where the influence of the relative crack length a/W on the distribution of the P/P 0 = f(v LL ) curves is illustrated, with appointed the yield point σ 0 and the strain hardening exponent n. Analyzing this graph, it can be argued with the fact that a specimen with a relative crack length a/W = 0.05 is characterized by the highest limit loadthe longer the crack, the lower the limit load P 0 . In Figure 14b, there is no reversal of the trend in the arrangement of the P/P 0 = f(v LL ) curves (with reference to the diagram presented in Figure 13b) the strain hardening exponent does not affect the limit load P 0 . Conversely, in Figure 14c, one can again observe a reversal of the tendency of the P/P 0 = f(v LL ) curves, which is due to the fact that the limit load strongly depends on the yield pointthe higher the yield point of the material, the greater the limit load level.
By transforming formula (2), it can be written that the η coefficient, numerically determined for energy A calculated from the area under the P = f(v LL ) curve, is equal to Based on formula (20), the values of the η coefficient were estimated for all 64 analyzed cases, and the results of the analysis were graphically presented in the form of graphs of changes in the η coefficient as a function of the J integral normalized by the product of the length of the noncracked specimen section b and the yield stress σ 0 ( Figure 15). The graphs present the changes in the η = f (J/(b·σ 0 )) curves and were prepared for three combinations to determine the influence of the relative crack length a/W (Figure 15a), the strain hardening exponent n in the R-O law (Figure 15b), and the yield point σ 0 (Figure 15c). As can be seen, the value of the η coefficient initially increases slightly with the increasing external load, and then, depending on the geometric and material configuration of the SEN(4PB) specimen, it tends to reach the saturation level or slightly decreases. It is noticeable that in the case of SEN(4PB) specimens bent with pure moment, in contrast to the SEN(B) specimens discussed in ref. [13], the value of the η coefficient strongly depends on the relative crack length a/W (Figure 15a), slightly from the strain hardening exponent n in the R-O law (Figure 15b), and is also insensitive to the change in the value of the yield stress σ 0 (Figure 15c). The value of the η coefficient clearly increases with the increase of the crack length a/W and very slightly with the decrease of the material hardening degree, expressed by the strain hardening exponent n in the R-O law. The value of the yield strength σ 0 does not affect the value of the η factor. These properties will be used in the next steps of the analysis of the obtained set of numerical results.
Due to the complex nature of the changes in the η = f(J/(b·σ 0 )) curves, it is not possible to mathematically describe the changes in the η coefficient as a function of the integral J normalized by the product of the yield point and the physical length of the noncracked section of the specimen (J/(b·σ 0 )). Hence, to simplify the analysis and provide a concise formula allowing to estimate the value of the η coefficient depending on the material characteristics and geometry of the specimen (expressed in this article by the relative crack length a/W), the analysis scheme shown in refs [37,38] was used, which was also used in ref. [13]. Let us normalize equation (2) by dividing it on both sides by the product of the noncracked section of the specimen b and the yield stress σ 0 . Such a procedure leads to the following dependence: Denoted by: The expression of the form is obtained which, after conversion, allows to write down the formula for the coefficient η: By plotting the changes of the normalized J integral denoted by J¯, as a function of the normalized energy A, denoted by Ā (Figure 16), and then analyzing them, we can confidently say that the value of the η coefficient is equal to the tangent of the φ angel (the slope of the ( ) = J f Ā¯curve to the abscissae). The analysis of Figure 7 leads to the same conclusions as in the case of the evaluation of the runs presented in Figure 15. The η coefficient depends on the relative crack length (a/W) and the strain hardening exponent (n). The value of the η coefficient is insensitive to the change in the yield stress of the material (σ 0 ). These relationships were used in the next stage of analysis to propose the simple formula to calculate the η coefficient depending on the relative crack length and the strain hardening exponent. Table 8 presents all numerical results obtained during analysis.
The value of the η coefficient does not depend on the yield point (slight changes in the value of the η coefficient, up to 5%, are observed for the extreme analyzed yield points σ 0 = 315 MPa and σ 0 = 1,500 MPa). As the crack length increases, the value of the η coefficient increases significantly. Conversely, a variable influence of the strain hardening exponent on the value of the η coefficient is observed. For n = 3.36 and n = 20 and in the case of very short cracks (a/W = 0.05), the value of the η coefficient differs by almost 100%. In the case of short cracks (a/W = 0.20), the difference is a maximum of 40%, and in the case of normative and very long cracks (a/W = 0.50 and a/W = 0.70, respectively), the difference is less than 6%. Figure 17a shows the total effect of the relative crack length a/W and the strain hardening exponent n in the R-O law on the value of the η coefficient, determined according to the diagram presented earlier (formulas (21)(22)(23)(24)(25)). As can be seen, the surface created on the basis of points determined by numerical calculations and analytical considerations proves the lack of dependence of the value of the η coefficient on the yield point and indicates a significant influence of the relative crack length and a weak influence of the strain hardening exponent in  the R-O law. Therefore, using the Table Curve 3D package [33], an equation approximating the surface built on the basis of 64 research points was adjusted, which takes the following form:   where A 1 -A 10 denote the coefficients of matching formula (26) to points on the η = f(n, a/W) surface presented in Figure 17a. The coefficient of determination R 2 for the presented figure is almost R 2 = 0.997, and the matching coefficients are presented in Table 9.
The analysis of Figure 17b shows that in selected cases, the approximation of all numerical results by the  formula (26) is not the correct solution. For selected geometric and material configurations, too large differences are obtained between the values of the η coefficient determined on the basis of numerical calculations and analytical approach, and the values obtained on the basis of formula (26). Therefore, to minimize the differences in critical points, it is recommended to use formula (27) where the A 1 -A 8 coefficients of approximation are presented in Table 10. The maximum difference in the case of using formula (27) and Table 10 is an error of about 7.5%.

Assessment of stress distribution in front of the crack tip
Full analysis and evaluation of stress distributions near the crack tip should include the tested 64 cases of material and geometric configuration, along with an analysis of the specimen reaching the limit load and full plasticization, which for elastic-plastic materials will not coincide with the limit load calculated according to the recommended formulas due to material strengthening. Such statements should be prepared for the successive main components of the stress tensor, which determine the level of known geometric constraints (Q stress, mean stress σ m , and triaxiality coefficient T z ). Due to too many For the purposes of this study, it was decided that the evaluation of stress distributions near the crack tip will be performed for selected geometric and material cases for the standardized external load P/P 0 ≈ 1.0. The following subjects were assessed: • effective stresses σ eff determined according to the Huber-Mises-Hencky hypothesis, normalized by the yield point σ 0 or by mean stresses σ m , calculated as xx yy zz m (28) where the component σ xx is the stress tensor component in the direction of the specimen thickness, and the components σ yy and σ zz are components in the crack propagation plane; • mean stresses σ m normalized by the yield point σ 0 ; • stress triaxiality factor T z , defined by Guo [21][22][23], calculated as yy zz (29) • principal components of the stress tensor -σ yy , σ zz , and σ xx .
The parameters mentioned earlier (effective stresses, mean stresses, and stress triaxiality coefficient) are considered by many researchers to be measures of stress triaxiality and geometric constraints, influencing fracture toughness and fracture development in elastic-plastic materials.

Figures 18 and 19
show the influence of the relative crack length on the mentioned parameters and selected main components of the stress tensor. The graphs were prepared for both the physical distance from the crack tip and the normalized distance, calculated as r·J/σ 0 . The farther from the crack tip, the lower the value of the stress tensor components (Figure 19), and thus the effective stresses or mean stresses. Along with moving away from the crack tip, the value of the triaxiality stress coefficient T z also decreases. As the crack length increases, the effective stresses and the stress triaxiality coefficient decrease are observed (Figure 18), with a simultaneous increase in the values of the main components of the stress tensor ( Figure 19)the shorter the crack, the lower the value of the successive components of the stress tensorit is visible at the analysis of graphs presented as a function of the normalized distance from the fracture tip. The same behavior is observed for the analysis of changes in the mean stress values σ m , also presented as a function of the normalized distance from the crack tip. It should also be noted that with the distance from the crack tip, in the tested measuring range, the value of the ratio of effective stresses and mean stresses increases -σ eff /σ m , while the increase in the crack length is accompanied by a decrease in its value (Figure 18a and b). Figures 20 and 21 show the effect of the hardening level on the stress distribution near the crack tip for SEN (4PB) specimens. The stronger the material strengthens (smaller value of the strain hardening exponent n in R-O law), the higher the stress tensor components are (Figure 21), similar to the effective stresses and average stresses normalized by the yield point ( Figure 20). However, as shown in Figure 20e  increase in the value of the stress triaxiality coefficient T z is observed and an increase in the value of the ratio of effective stresses and mean stresses -σ eff /σ m , which in many scientific papers is considered a natural measure of geometric constraints in fracture mechanics. The less the material strengthens, the greater the value of the quotient σ eff /σ m , as in the case of the triaxiality coefficient T z (Figure 20a and b). Conversely, Figures 22 and 23 show the influence of the yield stress on selected parameters considered to be measures of geometrical constraints in the fracture mechanics ( Figure 22) and on the main components of the stress tensor ( Figure 23). In the case of the main components of the stress tensor, the influence of the yield stress is insignificant when the stress distributions are considered as a function of the physical distance from the crack tip (Figure 23(a, c, and e)). However, the change of the physical distance into its normalized form, denoted by ψ = r·J/σ 0 , indicates that with the increase in the yield point, the value of the successive components of the stress tensor decreases (Figure 23(b, d, and f)). In general, the fastest changes are observed for the component opening the fracture surfaces, denoted by σ zz (Figure 23d). Also in the case of graphs presenting selected measures of geometric constraints (Figure 22(a and c)), a slight influence of the yield stress on the value of these parameters is observed (we are talking about the quotients σ eff /σ 0 , σ eff /σ m , σ m /σ 0 , and the T z parameter), as long as they are considered as a function of the physical distance from the crack tip, denoted by r.
Conversely, when considering the distributions of the quotients σ eff /σ 0 , σ eff /σ m , σ m /σ 0 , and the T z parameter as a function of the normalized position in front of the crack tip ( Figure 22(b, d, and f)), it should be noted that the SEN(4PB) plane strain W = 40mm B=1m E = 206GPa Q = 0.30 n = 10 a/W = 0.70  yield stress has a significant influence on the value of mean stress (Figure 22d) and the T z parameter (Figure 22f). The greater the yield point, the lower the value of the mean stresses σ m /σ 0 and the T z parameter. The value of effective stresses normalized by the yield point also slightly decreases (Figure 22b), and when normalizing the effective stresses by mean stresses, it can be noticed that the value of the quotient σ eff /σ m increases slightly with the increase of the yield stress of the material. The conclusions and observations presented in this section are generally reflected in the other cases that were considered during the preparation of this article. In the future, the author intends to extend the scope of the research work with the analysis conducted for cases of three-dimensional beams subjected to bending in the fracture plane with a pure bending moment.

The Q stresses as a measure of in-plane geometric constraints
The Q stresses in the form known and used today were defined by O'Dowd and Shih in the early nineties of the last century [24,25]. It was one of the attempts to improve the description of the stress fields in front of the crack tip in elastic-plastic materialsso that the theoretical distribution differed as little as possible from the real one and would introduce new elements to the fracture mechanics. It turned out that the introduction of a two-parameter description (enrichment of the HRR solution by taking into account the influence of in-plane constraints) improves the result; however, the notion of a plane strain state or a plane stress state is still used. turned out to be a simple approach that was used in solving engineering problems in the form of inclusion in European programs SINTAP [SINTAP, 1999] and FITNET [FITNET, 2006]. In a series of papers [24,25,41], the authors concluded that the results obtained with the FEM are accurate and compared the differences between it and the HRR field. They proposed a description of the stress field in the form of the formula: where the first part of the equation is the HRR solution [19,20], and the second part takes into account the influence of all other terms of the asymptotic expansion.
In equation (30), J is the J integral, σ 0 is the yield strength, ε 0 is the strain in the uniaxial tensile direction corresponding to the yield stress (ε 0 = σ 0 /E, where E is the Young's modulus), n and α are material constants from the Ramberg-Osgoog relationship, I n is a quantity dependent on the material through the strain hardening exponent n and the method of loading and thickness of the specimen, ( ) σ θ ñ , ij and ( ) σ θ n , ij are functions of the θ angle and the material constant n, which is similar to I n function determined numerically, q is the exponent of the second term of the asymptotic expansion, and Q results from the fitting of equation (30) to the numerical solution obtained using the finite element method (FEM). Conducting a series of studies, O'Dowd and Shih proposed the following equation to describe the stress field in front of the fracture front in elastic-plastic materials: which significantly simplifies the description of the stress fields in front of the crack tip in elastic-plastic materials. The undoubted advantages of the solution represented by formula (31) over many descriptions of mechanical fields using a multiterm expansion are its simple form and easy way to obtain the value of the Q parameter, which is often referred to in the literature as "Q stress." It should be remembered that in this case the Q stresses are not the second element of the asymptotic expansion, but the quantity taking into account influence the stress distribution of all terms of the higher-order expansion.
To avoid ambiguity in determining the value of the Q parameter, the authors in refs [24,25] determined that the most appropriate place to measure the Q stresses would be a point located at the distance r = 2.0·J/σ 0 in the direction of θ = 0. The choice of the direction in which the Q parameter value is measured is not accidental, but dictated primarily by practical reasons. O'Dowd and Shih [24,25] where (σ θθ ) FEM is the stress value determined numerically using FEM and (σ θθ ) HRR is the stress value resulting from the HRR solution.
The Q parameter also referred to in the literature as Q stress has many characteristics. In the case of the dominance of the plane stress state, the values of Q stresses are close to zero, while for the case of the plane strain state, the Q parameter assumes negative values [8,34,35,42,43]. Its value for the plane strain case domination depends on • geometry of a construction element (specimen)type of a structural element; • dimensions of a structural elementwidth and length of the crack; • materialyield point, strain hardening exponent; • method of loadingtension, bending; • the level of external load, expressed, for example, by the J integral (crack driving force) or by means of a real external load, normalized by a limit load. The Q parameter is quite commonly used in problems in the field of fracture mechanics. It is used not only to describe stress fields but also to assess the strength of structural elements containing defects using failure assessment diagram (FAD) diagrams to lower the level of conservatism [7]. The Q parameter is also used in determining the actual fracture toughness, if the appropriate fracture criteria are used for this purpose [33].
That is why it is necessary to know the Q stressesas a parameter correcting the stress field in front of the crack tip, as well as a parameter used in determining the actual fracture toughness. As a part of this study, the Q stress values were estimated for the considered beams with a crack subjected to four-point bending. This value was determined in the normalized distance (r·σ 0 )/J = 2, for the direction θ = 0, according to formula (32). The full results are 128 curves presenting the distribution of changes in the Q stress value as a function of the J integral (which is considered as the crack driving force) and the distribution of changes the Q stress as a normalized function the abscissa coordinate, written as log(J/(a·σ 0 )). Figures 24-26 present the selected results of the numerical analysis in a graphic form.
As can be seen, as the external load expressed by the J integral increases, the Q stress value decreases, assuming more and more negative values (Figures 24-26). The shorter the crack, the faster the changes in the Q stress value. When assessing the effect of the crack length, the behavior of the Q = f(J) and Q = f(log(J/(a·σ 0 ))) curves cannot be considered as fixedcharacter of these changes is different for specimens with short cracks (a/W = 0.05 and a/W = 0.20) and long cracks (a/W = 0.50, a/W = 0.70). The analysis of the impact of the crack length on the distribution of the Q = f(J) and Q = f(log(J/(a·σ 0 ))) curves is left to the reader based on Figure 24. It should be noted that the presented nature of the changes is generally consistent with the results obtained for other geometries with a predominance of bending, which the author mentioned in refs [8,35].
The analysis carried out in this article shows the unambiguous influence of the strain hardening exponent n ( Figure 25). As can be seen, the stronger the material strengthens, the lower the Q stress value for the same J integral level (Figure 25). The increase in the crack length is accompanied by faster changes in the Q = f(log(J/ (a·σ 0 ))) curvessee Figure 25c and e. Depending on the geometric and material configuration, the difference between the Q stress values for two extreme materials with different levels of hardening (n = 3.36 and n = 20), with the same J integral level, may be greater than 0.8.
The established dependence and the established influence on the distribution of the Q = f(J) and Q = f(log(J/ (a·σ 0 ))) curves are also observed when assessing the dependence of these curves on the yield point ( Figure 26). The higher the yield point, the higher the curves Q = f(J) and Q = f(log(J/(a·σ 0 ))) lie, and the Q stresses at the same level of the J integral assume higher values. As can be seen, the pattern of the Q = f(J) and Q = f(log(J/(a·σ 0 ))) curves changes with the increase of the crack length (Figures 25 and 26). The analysis of this fact, as in the case of Figure 24, is left to the reader. When assessing the value of Q stresses for two extreme yield stresses (σ 0 = 315 MPa and σ 0 = 1,500 MPa) considered in the course of numerical calculations, it can be noticed that at the same level of the J integral, the differences between them may be greater than 1.1. In the future, the author of this article intends to prepare a proprietary computer application that will allow viewing the collected Table 11: Coefficients of matching formula (33) to the numerical results for the estimation of the Q stress values, for SEN(4PB) specimens dominated by the plane strain state, with the relative crack length a/W = 0.05 results of numerical calculations in the graphic form and libraries in the form of PDF files. Selected combinations of Q = f(log(J/(a·σ 0 ))) curves presented in this article, along with other obtained results, were used to develop a catalog of numerical solutions, which were approximated by a third-degree polynomial. This approach allows to estimate the Q stress level without the need for tedious numerical calculations and quite complicated analyzes after their completion. The function selected to approximate the Q stresses is a function of the variable in the form of the parameter log(J/(a·σ 0 )): where A 1 , A 2 , A 3 , and A 4 are approximation factors, J is the J integral, a is the crack length, and σ 0 is the yield point. It is recommended to enter physical quantities into the aforementioned formula after they have been converted to values corresponding to the basic units. The use of formula (33) requires the user to know the current value of the J integral, the crack length a, the yield stress σ 0 , and the stain hardening exponent n in the R-O law. Formula (33) will correctly estimate the Q stresses in terms of material and geometric characteristics considered in this article. It should be remembered that the given values were obtained for a reference thickness of B = 1 m. The values of the approximation coefficients for the geometric and material configurations used in the research program, along with the coefficient of determination, are presented in Tables 11-14.  (33) to the numerical results for the estimation of the Q stress values, for SEN(4PB) specimens dominated by the plane strain state, with the relative crack length a/W = 0.20

Summary
This article presents the results of analytical and numerical research focused on the determination of selected fracture mechanics parameters for beams containing a crack in the state of four-point bending with the dominance of a plane strain state. Based on a wide program of numerical calculations carried out for 16 hypothetical elastic-plastic materials differing in the yield point and the strain hardening exponent in the Ramberg-Osgood law, details of numerical modeling using the finite element method are presented. In the further part of this article, devoted to the obtained results of numerical calculations, the influence of geometry and material characteristics on numerically determined values of the J integral, crack tip opening displacement, and load line displacement would be discussed. For these values, new hybrid solutions were proposed, allowing to estimate them without the necessity to carry out tedious numerical calculations.
In the next stage, the mutual relationship between the J integral and the crack tip opening displacement was discussed, taking into account the influence of the inplane constraints on the value of the d n coefficient defined by the Shih relationship, giving appropriate approximation formulas. This article also verified the Landes and Begley formula, originally used to determine the J integral based on the graph of force versus the load line displacement, also taking into account the influence of the inplane constraints on the value of the η coefficient, by making its value dependent on the yield point and relative crack length. Also in this part, the relevant numerical calculations were approximated, as presented earlier.
In addition to the assessment of the basic quantities in the field of fracture mechanics, this article presents a discussion on the stress distributions in front of the crack tip, as well as parameters commonly considered to be measurements of geometric constraints, among which a lot of attention was paid to the Q stresses defined by O'Dowd and Shih, with the impact of geometry and material characteristics on their values, as well as approximation formulas for 64 cases of geometric-material configuration, allowing to estimate their value without the need to carry out tedious numerical calculations.
The basic conclusions that can be drawn on the basis of the obtained results are as follows: • the lower the material hardening degree (higher value of the strain hardening exponent n in the R-O law), the higher the values of the elastic-plastic fracture mechanics parameters are obtained at the same level of external load; • the increase in the material strength (expressed as yield point σ 0 ) is also accompanied by an increase in the value of the J integral, crack tip opening displacement δ T , and load line displacement v LL at the same level of external load; • the values of the functions h h , 1 ⁎ 2 ⁎ and h 3 ⁎ (which are needed to determine selected fracture parameters using new hybrid solutions) depend on the strain hardening exponent n, the relative crack length a/W, and the external load, but it is not possible for the aforementioned three functions define equal dependence on the specimen geometry or material characteristics; • the significant influence of the strengthening exponent on the level of the h * function is observed;