*Article* **Change in the Properties of Rail Steels during Operation and Reutilization of Rails**

**Kassym Yelemessov 1, Dinara Baskanbayeva 1, Nikita V. Martyushev 2,\*, Vadim Y. Skeeba 3, Valeriy E. Gozbenko 4,5 and Antonina I. Karlina <sup>6</sup>**


**Abstract:** The paper considers the possibility of reusing previously used railway rails. The analysis is conducted using the standards and operating conditions of the rails of one of the Central Asian states, Kazakhstan, as an example. The operation of these rails causes significant stresses, while the surface layers are strengthened as a result of cold hammering. These phenomena significantly change the physical and mechanical characteristics of rails. As a result, they may not be suitable in terms of parameters for basic use but can be suitable for installation on other tracks. The conducted studies have shown that when the standard service life of the RP65 rail expires, the surface layer is deformed to a depth of up to 300 microns, hardness increases, and internal residual stresses are formed. These changes lead to an increase in the strength properties of the rails. However, at the same time, cracks originate in the surface layer of the rail, thus worsening operational characteristics. The RP65 rails are used under a cyclic load of 700 kN (which is determined by the national standard), withstanding 790,000 cycles. When the load is reduced to 510 kN, these rails can withstand the 2,000,000 cycles required by the standard without failure. Thus, these rails can be reutilized only on non-loaded and non-critical sections.

**Keywords:** rails; rail steel; reutilization; fatigue life; rail reliability

Reutilization of Rails. *Metals* **2023**, *13*, 1043. https://doi.org/10.3390/

Baskanbayeva, D.; Martyushev, N.V.; Skeeba, V.Y.; Gozbenko, V.E.; Karlina, A.I. Change in the Properties of Rail Steels during Operation and

Academic Editor: Haitao Cui

**Citation:** Yelemessov, K.;

Received: 20 April 2023 Revised: 20 May 2023 Accepted: 27 May 2023 Published: 30 May 2023

met13061043

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Rail transport plays a crucial role in the economies of the majority of countries, providing a significant share of all cargo turnovers and the transportation of most export and transit goods. This share in Central Asian countries may reach 50% or more. Railway rails perform an important function: the elastic transformation and transfer of loads proceeding from wheels to tracks, directing the running gears of locomotives and cars [1]. High requirements are imposed on the quality of rails during the operation of the railway track [2,3]. R65 rails are used in the manufacture of rail tracks in the most massive sections of traffic (main tracks) in Central Asia, such as the sections for cargo and passenger traffic. The speed of passenger trains in these sections reaches its maximum. In addition, freight trains are also involved in operations [4]. In such areas, steel with sufficiently good plastic properties and relatively low hardness (HB = 360–380 MPa) is used as the rail material. These rails are made of low-alloy carbon steel with a ferrite-pearlite structure [5]. There are many other sections in addition to the main railway, including industrial sections, safety dead ends,

crossovers, turnouts, marshalling, engine tracks, and station tracks, that have completely different requirements for rails [6,7]. Of the listed types of tracks, industrial ones make up a considerable proportion of their total number (ranking second as compared to the main tracks). Rails of the RP65 type are used on these tracks. Their difference from R65 rails consists of the type of material used in their production. The chemical composition of the RP65 material differs by having a higher content of sulfur and phosphorus as compared to that of R65. Owing to this fact, R65 has a higher plasticity and a lower tendency to crack. The same low-alloy carbon steel with a ferrite-pearlite structure is used to produce RP65 rails. The proportion of doped components in such steel is slightly lower, and the number of sulfur and phosphorus impurities is slightly higher, though such steel also has a slightly higher hardness and lower ductility compared to R65 steel [8]. The same material is often used for RP65 rails and R65 rails when the controlled parameters of the rail bed technological process decrease, thus changing the mechanical properties of the steel. In this case, stresses are formed in the surface layer during the long-term operation of the main track rail. The studies show that these are compressive stresses that will result in the rail surface hardening. The studies revealed that the stresses in each of the rail elements (head, wall, and base) are distributed individually [9]. These stresses depend not only on the operating conditions but also on the rail manufacturing technology [10].

Another phenomenon observed during the operation of rails that leads to a change in the properties of the surface layer is cold hammering. The authors of the work [11] show that in the subsurface layer of the rail head (located at a depth of 2–10 mm), the most significant physical strengthening mechanism is dislocation. This mechanism is due to the interaction of moving dislocations with stationary ones. In the surface layer of the rail head, however, the substructural mechanism is quite different. It is caused by the interaction between dislocations and small-angle boundaries of fragments and nanometer-range subgrains. The authors noted that tonnage ranging from 691.8 to 1411 million tons passes along the rails, increasing their strength 1.5–2.0 times. The authors of the work [12] also state that the main influence on the strengthening of rail steel is the dislocation substructure, which is formed during the rail operation. This work also notes that there is a change in the chemical composition of the rail surface during operation, which also contributes to the service life of the rails.

Both cold hammering and compression sets form in the surface layer of the rail and lead to its strengthening. As a result of these phenomena, after a certain period of rail operation, the strength properties of the rail surface will be higher than those in the initial state, i.e., before the operation. However, at the same time, cracks will form in the surface of the rail, leading to a decrease in its performance. The authors of the study [13] show that cracks often originate in the spots of various inclusions and defects in the surface layer of rails. Often, such cracks begin to develop on the side surface of the rail head. Grinding can prolong the life of the rails by crushing the head, removing damaged material, or moving the contact tape away from the stress concentrator [14,15]. There are very few research articles devoted to a comprehensive study of the properties of rails at various stages of their operation due to the significant labor input required. In addition, rail manufacturing standards and test equipment in different countries differ significantly from each other [16]. Although mass elements are similar, they may have different materials and slightly different geometric parameters. Therefore, the rails of each standard require a particular study. Regional operating conditions are also important. In [17], the authors study the influence of climatic changes on the track twist and rail joints under the conditions of the Spanish climate. The studies showed that temperature fluctuations strongly affected the durability of rails. Having analyzed the existing experimental data, the authors of [18] come to similar conclusions. In addition, based on the data obtained from the Swedish Portal for Climate Change Adaptation [19], they show that in climate conditions in some countries characterized by frequent variations in negative and positive temperatures during the year, the problem of the rail twist is critical.

Hence, it can be concluded that each region must conduct a specific study to accurately determine the changes in the properties of rails during operation. The results obtained through such studies allow us to think only about the possibility of reusing rails.

The aim of the work is to study comprehensively the properties of the high-duty rails of the R65 type after the warranty period of their operation and obtain information about changes that take place in their structure and properties. Obtaining comprehensive knowledge about such changes fills in the blanks in this knowledge area and helps determine the correct area for reusing these rails.

#### **2. Methods and Materials**

The R65 rails that were studied in this work had an actual operating time for the handled tonnage that amounted to 494 million t-km gross. The samples made from the RP65 rails manufactured by ARBZ LLP (Aktobe rail plant LLP, Aktobe, Republic of Kazakhstan) were used for comparison. The R65 rails were used on the main tracks. The choice of the RP65 rails for comparison is dictated by the fact that the rails are used to travel along industrial tracks for less critical purposes. The R65 rails with a hardened surface will be considered as an option to replace the RP65 rails that are less crucial.

The chemical composition of the materials used for the rails is shown in Table 1. Table 1 reflects the chemical composition of the used rails according to the industry standard. The chemical composition of the R65- and RP65-grade rail steels is almost identical, with the content of alloying elements (Mn, Si, and V) differing by no more than 0.1%. Such a difference does not significantly and noticeably influence the mechanical properties. The main difference is observed in the content of harmful impurities such as sulfur and phosphorus. Even minor changes in the concentration of these substances can significantly worsen the plastic properties of the steel.


**Table 1.** Chemical composition of the rail steel.

The metallographic specimens were made from the samples obtained using conventional technology to perform structural studies. To identify the carbon steel structure, a 3% solution of nitric acid in ethyl alcohol was used. The surfacing layer structure was detected by means of a solution of HNO3 and HCl acids in a ratio of 1:3. Structural studies were carried out using a Carl Zeiss AxioObserver Z1m light microscope and a Carl Zeiss EVO 50 XVP scanning electron microscope (Jena, Germany). The phase composition was studied using an ARL X'TRA (Thermo Fisher Scientific, Waltham, MA, USA) X-ray diffractometer in the CuKα radiation.

The following methods were used to identify defects in the welded joints: a visualoptical method using a Carl Zeiss AxioObserver A1m (Carl Zeiss Microscopy Deutschland GmbH, Oberkochen, Germany) microscope; a capillary method; and an eddy current method using a VD–70 eddy current flaw detector.

#### *2.1. Mechanical Testing of Samples*

The hardness of the rolling head surface was determined as follows: The rail's hardness was monitored using a Brinell device. The hardness of the surface and cross section of the rail was determined by a WILSON BH3000 (ITW Test & Measurement GmbH, Duesseldorf, Germany) hardness tester according to GOST 9012-59 (Figure 1). The spot for determining the hardness of the rail's rolling surface was cleaned to remove scale and a decarbonized

metal layer to a depth of no more than 0.5 mm. The roughness of the cleaned surface was no more than 25 μm.

**Figure 1.** Measuring the hardness of the samples with the WILSON BH3000 hardness tester.

The mechanical tension properties of the rail were determined using the Instron HDX-1000 universal rapture test machine. The blanks of samples for the tensile test were cut along the rolling direction taken from the upper part of the head in the fillet zone as close as possible to the surface at a distance of at least 150 mm from the end of the rail (Figure 2). The main dimensions of the samples are also shown in Figure 1. The mechanical characteristics were recorded automatically and calculated by the control and calculation software of the test machine. The measurement error of this machine is ±0.5%. Tensile strength mechanical characteristics were determined according to the uniaxial tension scheme on the cylindrical samples (Figure 2), which had a pick-up movement rate of 10 mm/min.

**Figure 2.** The spot where the samples for tensile testing were cut and the sample dimensions. (**a**)—the place where the sample was cut from the rail, (**b**)—test specimen drawing.

The experimental results were statistically processed in Statistica (StatSoft Inc., Tulsa, OK, USA), Table Curve 2D, and Table Curve 3D software. For each value (each point in the diagram), the tests were conducted on at least 5 samples.

#### *2.2. Cyclic Life*

The cyclic life during fatigue tests is determined according to GOST 25.502 under hard loading (strain control) of the samples with a constant amplitude of the full strain (longitudinal) equal to 0.00135. Cyclic durability tests were performed using the Instron 8801 (Instron Engineering Corporation, Norwood, MA, USA) universal machine (Figure 3). This machine has built-in sensors for measuring the deformation value. Sample blanks for the test were cut along the rolling direction from the upper part of the head in the fillet area as close as possible to the surface at a distance of at least 150 mm from the end

of the rail (Figure 2). The basic dimensions of the samples are also shown in Figure 1 (Loading scheme—cyclic tension-compression). Cyclic durability tests were conducted indoors at an ambient temperature of 22 ◦C and a relative humidity of 65%. The samples had a temperature of 22 ◦C before the tests. A longitudinal uniaxial cyclic load was applied to the sample, which had a loading cycle asymmetry coefficient of minus 1 and a loading frequency of 40 Hz. The test base is 5 million loading cycles. The testing was terminated when a crack or fracture in the sample formed or when the test base was reached.

**Figure 3.** Appearance of the test complex, Instron 8801. 1—hydroelectric power station, 2—controller, 3—computer, 4—power frame, 5—control panel, 6—columns, 7—rigid horizontal traverse, 8—force action sensor, 9—sample, 10—hydraulic grips, 11—hydraulic servo drive column, 12—panel for controlling the hydroelectric station, and 13—valves for controlling the stiffening beam position.

#### *2.3. Rail Endurance Limit (Fatigue Tests)*

The tests were carried out on full-profile samples (1200 + 10 mm long), cut from the rails by cold mechanical cutting methods. The loading scheme is a flat three-point symmetric bending (Figure 4). The distance between the lower supports is 1000 ± 5 mm. The upper punch is installed in the middle between the supports within 500 ± 5 mm. The samples are tested under soft loading (force control) in the position "head down" of the rail when the asymmetry of the loading cycle is plus 0.1. The test base is 2 million cycles.

**Figure 4.** Experimental plant and a general view of tests for determining the endurance limit of the rails.

Cyclic life tests are performed indoors at an air temperature of 22 ◦C and a relative air humidity of 65%. The temperature of the samples before testing was 22 ◦C.

Measuring tools or equipment must provide a cyclic load of at least 1000 kN and have a maximum permissible relative measurement error of ±3%. The maximum relative error of the sensor must be ±2%. The loading frequency is 10 Hz, with a maximum relative error of ±2%.

#### *2.4. Residual Stresses*

Residual stresses in the rail neck are determined by the divergence of the groove as a rail height difference (H2-H1) along the axis at the end of the sample before and after cutting the groove according to the scheme shown in Figure 5. The samples intended for studying residual stresses by the mechanical method were prepared on an electrochemical die-sinking and hole-contouring machine, SFE-12000M (STANKOFINEXPO, Kirov, Russia). Such equipment was chosen because no additional technological residual stresses were introduced when using this method of preparing experimental samples. Initially, a piece of the rail having a length of 600 ± 3 mm was cut off, and then a groove 6 ± 1-mm wide was cut through, and the accuracy of this size was ensured by the accuracy of manufacturing with an electrode tool.

**Figure 5.** Schematic of residual stress control in a track neck. (**a**)—rail before testing, (**b**)—rail after cutting for residual stresses testing.

The control of residual stresses in the rail bottom is referred to as "periodic tests" and, according to item 6.5 "Periodic tests" of GOST 51685-2013, is carried out at least once every three years. The level of residual stress in the rail bottom is determined from the fullprofile samples. The samples, 1.0 ± 0.1-m long, were cut from six rails by cold mechanical cutting. It is necessary to perform abrading to a depth of 0.3 to 0.5 mm on the bearing surface of the rail bottom in the middle part of the sample, after which a strain gauge is attached to the cleaned area in the longitudinal direction according to the recommendations of the manufacturer of the sensor. The permissible relative error of the sensor shall not exceed ±1%.

#### **3. Results and Discussion**

#### *3.1. Study of Rail Microstructure*

First, the microstructures of already-used and new rails were studied. The photos of the etched microstructure of the used R65 rails and the new R65 rails are shown in Figure 6.

**Figure 6.** Microstructure of the studied R65 rails: (**a**) surface of the rail after the standard service life of 494 million t-km gross; (**b**) new (unused before) rail in a state of delivery; (**c**) pearlite grain on the rail surface after the standard service life; and (**d**) ferrite-pearlite structure section of the rail surface after the standard service life (at a depth of ~3 mm from the surface).

The metal microstructure of the rails after operation is represented by a lamellar pearlite estimated at 2-3 points (0.3–0.4 microns) with scattered ferrite sections along the grain boundaries. The number of ferrite grains present in the used R65 rails (Figure 6d) is less than 5%, and the average size is about 0.24–0.26 μm. There is practically no perlite in the new R65 rails (Figure 6b), and bainite was not detected in the microstructure of the studied rail. When moving away from the surface of the rail rolling head, perlite acquires a more coarse-grained structure. This is related to the specificity of rail manufacturing technology. The rails taken for this study were made by hot rolling. Since the rail cools unevenly after rolling, a finer-grained structure is formed on the surface; however, when moving away from the surface of the rail head, the pearlite grain size becomes larger. Characteristic images of this structure with a larger grain are shown in Figure 6d.

A deformed structure is observed to a depth of up to 300 microns on the etched cuts from the surface of the working fillet. The amount of the decarbonized layer on the surface identified by the solid ferrite grid does not exceed 320 μm (Table 2). Local light-etched areas of a work-hardened metal are observed on the rolling surface, along which cracking develops (Figure 6a).


**Table 2.** Analysis of microstructure test results.

According to the classical definition, pearlite is a structural component of steel and cast iron, a eutectoid mixture of ferrite and cementite [20]. When crystallized under normal conditions, pearlite has a lamellar structure consisting of alternating plates of pearlite and cementite [21]. The studies performed in the present work revealed numerous imperfections in the structure of lamellar pearlite, namely, an alternating "comb" type structure (Figure 7a) and breaks of cementite plates (ferrite bridges) (Figure 7b). Quite often, there are curved plates of cementite of varying thickness in the pearlite colony.

**Figure 7.** Pearlite structure in the surface layer: (**a**) alternating pearlite structure and (**b**) broken structure of cementite plates.

The ferrite plates of pearlite colonies quite often have an alternating light and gray contrast. This diffraction contrast indicates that the ferrite plates are obviously divided into weakly misaligned areas as a result of elastic stresses.

The metallography of etched cuttings (Figure 6d) revealed that structurally free ferrite grains are present in the steel structure of the used rails, i.e., ferrite grains, in the volume of which there are no carbide phase particles. The relative content of such grains is small and does not exceed 5% of the steel structure. Ferrite grains, as a rule, are located in the form of interlayers along the borders of pearlite grains (Figure 6a) or at the joints of the borders of pearlite grains (Figure 6d).

Since the rails are obtained by means of rolling, the absence of line defects is important. Therefore, we conducted a rail study for the parameters of individual globular defects (ED) and line defects (EB). To identify the possibility of reutilizing rails, it is necessary to study their structure for the presence of non-metallic inclusions since such inclusions weaken the properties of steel. In general, only the shape and size of inclusions may change during operation, and the number of inclusions may only vary in the surface layer. To determine the average sizes of the inclusions, 8 samples for the new RP 65 rail and the already exploited R65 rail were studied. Eight samples were taken from different rails of each type to reproduce the results. A non-etched microstructure study showed that the average diameter of the inclusions for the new and already used rails is almost the same (Table 3). However, the length of inclusions for the latter is much greater, and at the same time, the total coefficient Ka for the used rails is still lower than the standard value.

#### *3.2. Study of the Mechanical Properties of Rails*

The studies of the mechanical properties of the rails showed that there was a slight difference in the values of the strength limits and the plastic characteristics of the new RP65 rails and the previously exploited R65 rails (Table 4, Figure 8). The obtained results demonstrated that the strength properties (σB) for the material of the previously used R65 rails are 40 MPa lower than those of the material of the new RP65 rails. However, at the same time, the plastic properties of R65 are higher by about 8.5%. This can be explained by an appropriate combination of factors. On the one hand, the previously used rails have larger stitch defects, which should lead to a decrease in plastic properties. On the other hand, the initial chemical composition of R65 contains less sulfur and phosphorus. In addition, the surface layer undergoes plastic deformation during operation and has a finer-grained structure.


**Table 3.** Analysis of steel contamination with non-metallic inclusions.

\* Note: non-metallic inclusions: ED—individual globular inclusions; EB—line globular defects.

**Table 4.** Tensile mechanical properties.


**Figure 8.** Tensile diagram for the studied samples cut from R65 and RP65 rails.

The difference in the hardness of the rolling surface along the length of the rails and samples was determined by three measurements on the middle line of the rolling surface. An interval of at least 25 mm was taken for each of the three samples, which were taken from the ends and the middle part of the rail or on the rail surface (Tables 5 and 6).

**Table 5.** Hardness of the head rolling surface.



The diagram of hardness measurements along the rail section is shown in Figure 9.

**Figure 9.** Hardness measurement: (**a**) schematic of hardness measurement and (**b**) dependence of hardness on the distance from the rail surface.

Long-term operation of rails is usually accompanied by deformation and transformation of the material structure [22]. This can be justified by the hardness values obtained for the operating rails. The upper hardness limit is HB401So for the R65 rails. The hardness at a distance of 10 and 22 mm from the rolling surface must be 345–360 (HB10 mm) and 325–350 (HB25 mm). The hardness of the studied rails is within the upper boundary of these values or exceeds them, and the property changes are consistent with the changes in the rail structure. The quantitative analysis of the morphological state of the steel structure performed in this work showed that the operation of rails is accompanied by the transformation of the state of lamellar pearlite grains, namely, the destruction of cementite plates. According to the photographs in Figure 6, regardless of the position of the material's analyzed volume (rolling surface or fillet surface), the destruction of the lamellar pearlite structure is at its maximum in the surface layer of the rails, with a thickness of less than 2 mm. However, the destruction degree of the lamellar pearlite structure depends substantially on the position of the volume to be analyzed; namely, on the rolling surface, the relative content of the broken pearlite grains is more than 2 times higher than the content in the surface layer of the working fillet.

As noted above, pearlite grains become broken during rail operation. One of the main mechanisms of such destruction during plastic deformation of the steel is cutting cementite plates by gliding dislocations [23]. The operation of rails is accompanied by an increase in the level of elastoplastic stresses in the steel. The value of elastic plastic stresses in the steel, in accordance with [24], is characterized by an excessive dislocation density and curvature-torsion amplitude of the material crystal lattice. Both of these characteristics of the steel are noted to be determined in the analysis of the bending extinction contours of the material.

The above results of the study of the rail metal after long-term operation indicate the structural transformation of the lamellar pearlite.

The change in the cementite's elemental composition during crushing is minimal. At the initial stage of the transformation, the cementite plates of the pearlite colony are covered with gliding dislocations. This is accompanied by breaking the cementite plates into separate, weakly oriented fragments. Then, the structure of carbide changes as the plastic deformation degree of the material increases due to the stripping of carbon atoms from the cementite crystal lattice. It is worth remembering that this process is possible due to a noticeable difference in the average bond energy of carbon atoms with dislocations (0.6 eV) and iron atoms in the cementite crystal lattice (0.4 eV) [25].

The considered deformation transformations of the rail steel structure during operation on the railway should be noted as not adversely influencing the product's cyclic life. Table 7 shows the results of testing the cyclic life of the rails. The tests were carried out on the rolling surface and, for comparison, on the cross section at a distance of at least 10 mm from the rolling surface. The tests were stopped when the sample became cracked or fractured or when the test base was reached.

**Table 7.** Analysis of the cyclic life test.


The test results are considered positive if there are no fractures or cracks in all the tested samples upon reaching the test base. The test results are considered negative if the formation of a crack or fracture in at least one sample occurred within a number of loading cycles less than the test base.

The studies of the hardness along the cross section of the used rail show that the changes in properties can be considerable. Hence, the rolling of the wheel on the surface of the rail not only causes changes in hardness but can also form residual stresses at a considerable depth. The testing of the samples to determine residual stresses showed that there are considerable differences in values for the new and already used rails (Table 8). Residual stresses in the neck of the used rails are 1.42 times higher than those in the new rails (Figure 10). However, the stress values are in the upper range of values allowed for 2 mm rail use. A similar situation is observed for the middle part of the rail base. The residual stress level for the used rails is 1.25 times higher than that for the new ones (Table 9). In addition, the residual stresses of the used rails lie in the upper range of permitted values for use. Such differences in voltage values are certainly related to the operating conditions of the used R65 rails. Rolling the wheels of railway trains under a significant load deforms not only the surface layer but also the rail. This, in turn, determines the appropriate level and nature of the distribution of residual stresses. In the new rails, small residual stresses are formed at the manufacturing stage. Since the rails are made by hot rolling, internal stresses are formed in them [22].


**Table 8.** Residual stress indicators in a rail neck.

**Figure 10.** Samples after residual stress control in a rail neck.

**Table 9.** Analysis of test results for determining residual stresses in the middle part of the rail base.


The tests show that during operation, there are changes in the structure of the surface layer. Hardness increases along almost the entire rail section, and residual stresses are formed, but the defective surface layer of the rail also increases. The number and size of individual defects do not change, but the length of horizontally extended defects increases. In addition, cracks will form on the surface of the rail head. In general, such changes significantly affect one of the main operational properties of rails, i.e., the endurance limit.

The normalized values for using the rails in Central Asia are 2,000,000 cycles under a load of 700 kN for rails R65 and RP65. The new rails certainly meet these requirements (Table 10). In the case of the previously used rails, the maximum load under which they could operate for 2,000,000 cycles was 510 kN. Under a standard load of 700 kN, the used rails can only withstand 790,000 cycles (Figure 11a).

**Table 10.** Analysis of test results for the determination of the rail endurance limit.


**Figure 11.** Rail failure after cyclic tests: (**a**) under a load of 700 kN, number of cycles: 790,056; (**b**) fracture under a load of 600 kN, number of cycles: 1,111,966.

This result is largely due to the formation of cracks in the surface layer of the rails. A photograph of the fracture surface (Figure 11a) shows that the fracture originated on the side surface of the rail. The reason for the crack development was the zone of lateral collapse of the rail head formed under operating conditions (Figure 12a). This crack nucleus led to its fatigue development. Figure 12a clearly shows the zone of stable fracture development. It is characterized by fatigue lines representing approximate concentric contours. The focus of these contours is at the origin of the fatigue crack (Figure 12b). The surface of this zone is smooth and level. This effect is the result of cold hammering caused by repeated presses of two surfaces of the crack on each other.

**Figure 12.** Failure of a used rail after endurance limit tests. (**a**) A general view of the rail head and (**b**) the place where the crack originated on the rail fluting.

A set of experimental studies showed considerable changes in the structure of the R65 rail during its operation. The hardness of the surface layer and the rail itself increases to a sufficiently significant depth (up to 35 mm). In addition, considerable stresses are formed in the rail. It should be noted that both hardness and stress increase to the upper limit of the values allowed by the state standard (sometimes slightly exceeding it). Despite the increase in hardness and the formation of stresses, the plastic and strength properties of the rails remain almost at the same level and do not fall outside the tolerance field of the values determined by the state standard. The change in hardness is conditioned by changes in the structure of the rail. When the structure is deformed, the condition of the grains of lamellar pearlite is transformed, with the maximum amount of structure damage occurring in the surface layer of rails with a thickness of less than 2 mm. The lengths of the line defects also change, which then affects the cyclic life of the rails. New RP65 rails under a standard load of 700 kN can withstand 2,000,000 cycles. The previously used R65 rails under the same load can withstand significantly less than 790,000 cycles. They can withstand the 2,000,000 cycles determined by the state standard only under a load of 510 kN.

The changes that took place in the mechanical properties, especially in the cyclic properties, are related to both the operation intensity, structure, and mechanical composition of the R65 and RP65 rails. The analysis of the chemical composition of the analyzed steels is shown in Table 1. It shows that the chemical composition of the main alloying elements differs slightly, by no more than 0.1%. Such changes in the composition vary the properties of the steel by no more than a percentage [26]. Moreover, manganese in the RP65 steel may increase by 0.1%, and vanadium may increase by less than 0.05%. Such changes in the chemical composition will not noticeably change the mechanical properties [27,28]. The main difference between these steels is the content of sulfur and phosphorus. In RP65 steel, the content of sulfur and phosphorus can be higher by 0.01–0.02%, which somewhat reduces its plastic properties [27]. Based on this particular chemical composition, the R65 rails are designed for the main, working tracks [28], while the RP65 rails are intended for industrial routes, where the intensity of the movement is lower and there is not as much frequent deforming effect. The initial mechanical properties of the new R65 rail are slightly higher than those of the R65 rail. However, our tests have shown that both brands can withstand a certain standard load. Our studies have shown that during operation, microcracks form in the surface layer of the R65 rail, and the surface layer is decarbonized. These processes run in parallel with the processes of riveting the surface layer and the formation of internal stresses in the rails. The result is that the mechanical properties of the R65 rail (hardness and tensile strength) practically do not change or even increase by a small amount (5–10%). The formation of elongated defects in the surface layer and the formation of microcracks during operation lead to a decrease in the cyclic properties (rail endurance limit) of the rail. As a result, while the main mechanical properties of the rail material remain at the required level, the value of the rail endurance limit indicator for the used R65 rails becomes lower than that of the RP65 and no longer meets the requirements of the standard.

Based on this fact, the R65 rails can be said to have been industrially used for a standard period, which is not suitable for reutilization in the Republic of Kazakhstan. Their use is possible only in low-congested areas or areas without intensive movement (storage areas, dead ends, and other inactive routes). The number of passageways per day determines the possibility of using these rails. The regulatory framework for testing, amounting to 2,000,000 cycles that the rail must withstand, is determined by a guaranteed 20-year period of operation [29]. At the same time, they assume that up to 90 trains should pass per day. Based on the obtained results, the maximum number of cycles that the previously used R65 rail can withstand is 790 thousand, which is 2.53 times less. Therefore, it is possible to use these rails only within the sections where the traffic intensity is lower than 35 trains per day on average or various sections with non-intensive traffic, i.e., stops, dead ends, and other inactive ways. At the same time, the already-used R65 rails will provide the operation for the required period of 20 years.

At the same time, the standards and requirements for rails employed in different countries may differ significantly, as determined by rolling stock weight standards, climatic conditions, etc. Standards and requirements for the number of cycles of suitable rails can be different. Therefore, the suitability of rails for different countries can vary, necessitating meeting the standards of the country to determine the possibility of rail reutilization.

#### **4. Conclusions**

1. The studies revealed that the microstructure of the rail R65 metal after operation is represented by a 2–3-point lamellar pearlite (0.3–0.4 μm) with isolated ferrite sections along the grain boundaries, while the number of ferrite grains in the previously used rails is less than 5% (average size is 0.24–0.26 μm). The depth of the deformed layer is up to 300 μm, and the value of the decarbonized layer of the surface, identified by a solid ferrite network, does not exceed 320 μm. The average diameter of non-metallic inclusions for the new RP65 and already-used R65 rails is almost the same. However, the length of inclusions for R65 is much higher. After the warranty period of R65

operation, the length of line-extended defects increases. The long-term operation of the rails is accompanied by deformation and transformation of the material structure. The hardness values of the previously used rail during the warranty period were at the level of HB 350–360.


**Author Contributions:** Conceptualization, K.Y. and D.B.; methodology, D.B.; validation, V.E.G. and A.I.K.; formal analysis, V.E.G. and A.I.K.; investigation, V.E.G.; data curation, K.Y. and D.B.; writing—original draft preparation, K.Y., N.V.M. and V.Y.S.; writing—review and editing, N.V.M. and V.Y.S.; supervision, N.V.M. and V.Y.S.; project administration, N.V.M. and V.Y.S.; funding acquisition, A.I.K. and V.E.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the Ministry of Science and Higher Education of the Republic of Kazakhstan within the framework of the BR18574141 project Theme of the project: "Comprehensive multi-purpose program for improving energy efficiency and resource saving in the energy and mechanical engineering for the industry of Kazakhstan", for which the authors express their deep gratitude to them.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available from the corresponding authors upon reasonable request.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Changlong Wen 1,\*, Yanbing Zheng 1, Dong Mi 1,\*, Zhengming Qian <sup>1</sup> and Hongjian Zhang <sup>2</sup>**

<sup>1</sup> AECC Hunan Aviation Powerplant Research Institute, Zhuzhou 412002, China


**Abstract:** Two shape optimization methods, based on non-parametric and geometric parameters, were developed to address stress concentrations in the vent holes of gas turbine flow guide disks. The design optimization focused on reducing the maximum equivalent stress at the hole edge in an aero-engine gas turbine flow guide disk. The effectiveness of both methods in achieving this objective was studied. The results indicated that the non-parametric-based optimization method reduced the maximum equivalent stress at the hole edge by 24.5% compared to the initial design, while the geometric parameter-based optimization method achieved a reduction of 20.2%. Both shape optimization methods proved effective in reducing stress concentrations and improving fatigue life. However, the non-parametric shape optimization method resulted in a better design for the vent holes based on the study's findings.

**Keywords:** shape optimization; flow guide disk; vent hole; design optimization

#### **1. Introduction**

Aero-engine components feature a significant number of hole features that serve various purposes, including connection, weight reduction, air system functionality, and heat transfer. These holes play a crucial role in the overall performance and efficiency of the engine. However, it is essential to carefully consider the location and design of these holes to prevent the occurrence of severe stress concentrations at their edges. Stress concentrations can lead to structural fatigue and a reduction in a component's overall lifespan [1–4], and the stress concentration problem has received great attention from aero-engine strength researchers [5,6]. Among the various components in an aero-engine, the gas turbine flow guide disk vent serves as a prime example, highlighting the challenges faced in the design of such component structures. The design of these structures is often constrained by the relative positions of the holes and the required cross-sectional areas. As a result, the design space for these components becomes limited, posing significant difficulties for engineers and designers.

To address these challenges and optimize the design of aero-engine components with hole features, researchers have undertaken extensive studies in the field of design optimization. These studies have aimed to enhance the performance and longevity of these components by reducing stress concentrations and improving fatigue life. Chen et al. [7,8] focused their research on a specific component—the high-pressure turbine disk-mounting side bolt hole. They approached the optimization problem by describing the hole boundary using a biaxially symmetric-shaped hole. By modeling and optimizing the geometric dimension parameters of this shaped hole, they successfully reduced the maximum equivalent stress value at the hole edge, thereby enhancing the component's performance and fatigue life, and they also expanded the biaxial symmetry non-circular hole to the uniaxial symmetry non-circular hole and obtained different results. Sun et al. [9] pursued a different approach by conducting shape optimization based on the hyper elliptic equation

**Citation:** Wen, C.; Zheng, Y.; Mi, D.; Qian, Z.; Zhang, H. Design for the Vent Holes of Gas Turbine Flow Guide Disks Based on the Shape Optimization Method. *Metals* **2023**, *13*, 1151. https://doi.org/10.3390/ met13071151

Academic Editor: Alireza Akhavan-Safar

Received: 29 May 2023 Revised: 19 June 2023 Accepted: 19 June 2023 Published: 21 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and the sequential response surface method for the open hole characteristics of a rotary casing. Their research demonstrated that shape optimization techniques could effectively improve the dynamic characteristics of the structure, leading to further enhancements in performance and reliability. In another study, Yan et al. [10] employed a multi-island genetic algorithm and introduced a non-circular-shaped hole design to optimize the gas turbine flow guide disk vent hole. By combining this approach with finite element analysis sub-model technology, they successfully reduced the maximum equivalent stress value at the hole edge and improved the component's fatigue life. In another article by Yan et al. [11], surrogate-based optimization with improved support vector regression was used for the non-circular vent hole on an aero-engine turbine disk. Hossein et al. [12–14] discussed the stress concentration factors (SFCs) in circular hollow-section X-connections retrofitted with a fiber-reinforced polymer under different load conditions, and the SCFs were greatly reduced. Zhu et al. [15] focused their efforts on optimizing the shape of a bolt hole in an engine sealing disk using multi-segment curves. The results of their research were impressive as they achieved a remarkable 32.8% reduction in the maximum equivalent stress at the hole edge compared to the original single round hole design. These optimization strategies have exhibited significant potential for enhancing a component's performance and durability.

In the current research landscape, the optimization of aero-engine component designs with hole features has gained significant attention. Researchers worldwide are striving to improve the structural integrity and performance of these components through innovative approaches and advanced optimization techniques. One prominent area of research has focused on advanced optimization algorithms and methodologies. Genetic algorithms, particle swarm optimization, response surface methods, and surrogate modeling techniques are some of the widely employed optimization strategies [16–21]. These algorithms have aimed to explore the design space efficiently and identify optimal solutions that reduce stress concentrations and enhance fatigue life.

Researchers have primarily relied on traditional geometric parameter-based optimization methods and finite element analysis to enhance the design of aero-engine components featuring hole features. The studies mentioned above have demonstrated the effectiveness of various optimization techniques in reducing stress concentrations and improving the overall performance and fatigue life of these structures. However, certain limitations, such as vent area optimization, still require further exploration. By addressing these challenges, researchers can continue to advance the field and contribute to the development of more efficient and reliable aero-engine components. The current approach employed by researchers for optimizing the design of hole structural features involves using the traditional geometric parameter-based shape optimization method. This method relies heavily on the engineering design experience of researchers and utilizes specific geometric parameters such as elliptical, multi-circular arc, or shaped holes to define the design boundaries. The optimization outcomes are highly dependent on the selection of these design variables. However, an alternative approach known as non-parametric shape optimization offers a promising solution [22–27]. Unlike the traditional method, non-parametric shape optimization does not rely on geometric curves to describe the design boundaries. This method holds the potential to overcome the limitations of poor design optimization results caused by the lack of relevant structural design experience among designers. It also offers a theoretical avenue to enhance the design optimization process.

To explore the effectiveness of the non-parametric shape optimization method, this study focused on the design optimization of a gas turbine flow guide disk vent hole. By utilizing the non-parametric approach, we aimed to obtain optimal design solutions that could outperform those achieved through the traditional geometric parameter-based shape optimization method. Furthermore, a comparative study was conducted between these two methods, providing valuable insights and serving as a reference for future design optimization endeavors involving similar structural features. By adopting the nonparametric shape optimization method, researchers can aspire to enhance the efficiency

and accuracy of design optimization for hole structural features, ultimately contributing to advancements in the field of engineering design and promoting innovative solutions for various applications.

#### **2. Shape Optimization Method and Process**

#### *2.1. Geometric Parameter-Based Shape Optimization Method and Process*

To ensure minimal and independent design variables, as well as facilitate ease of processing and manufacturing, researchers employing the geometrical parameters of the shape optimization method must carefully choose a suitable curve to construct the parametric model for the hole boundary. Common approaches for constructing the hole scheme include selecting straight lines, multi-segment arcs, and ellipses.

In this paper, we present the flowchart of the geometric parameter-based shape optimization method (illustrated in Figure 1). Initially, a CAD parametric model is established, where the geometric parameters are modified and linked to the CAD model. This model undergoes preprocessing through finite element software, and it is subsequently solved. The optimizer then verifies the constraints and determines if the optimization results are optimal. If the design fails to meet convergence criteria, the process returns to updating the design by modifying the CAD geometric parameter model, and the aforementioned steps are repeated. It is important to note that this optimization method requires re-meshing after each iteration, and the selection of design variables significantly impacts the resulting optimization outcomes.

**Figure 1.** Flowchart of geometric parameter-based shape optimization.

#### *2.2. Non-Parametric-Based Shape Optimization Method and Process*

Non-parametric shape optimization is a method of optimizing designs that does not rely on specific geometric parameters to describe the shape. Instead, it offers flexibility by allowing the shape to be modified freely within the design space. This method is particularly useful in the finite element analysis (FEA) framework, where mesh deformation techniques can be employed to alter the mesh geometry.

In non-parametric shape optimization, the shape is typically represented using a set of control points or a mesh. The position of these control points or nodes can be adjusted to deform the shape. Mesh deformation techniques, such as the Free-Form Deformation (FFD) or Radial Basis Function (RBF) [28,29], can be utilized to smoothly modify the shape while maintaining mesh quality.

The optimization process involves defining an objective function that quantifies the desired performance of the design. This function can be based on various criteria, such as minimizing stress concentration, maximizing stiffness, or optimizing fluid flow characteristics. Additionally, constraints may be included to satisfy design requirements or limitations.

To optimize the shape, mathematical optimization algorithms, such as gradient-based or evolutionary algorithms, are employed. These algorithms iteratively adjust the positions of the control points or nodes to improve the objective function while satisfying the constraints. The optimization process continues until a satisfactory design solution is obtained [22].

Using the non-parametric shape optimization method, we can consider a hole-type structure as an example. In this approach, the coordinates of finite element nodes on the hole boundary are selected as design variables for optimization. During the optimization process, each node can be individually moved from its adjacent nodes. A schematic representation of coordinate movements for the design variable nodes is shown in Figure 2. Additionally, the nodes near the design variables can be adaptively moved, and to prevent jagged boundaries in the optimized structure, a mesh-smoothing algorithm is employed to ensure mesh quality.

**Figure 2.** Schematic of design node movement for a hole-type structure.

In this study, we illustrated the non-parametric-based shape optimization process, as depicted in Figure 3. A finite element model of a vent hole was set up using the commercial FEA tool, ANSYS 19.2(ANSYS Inc., Pittsburgh, PA, USA). The non-parametric optimization process was performed by Tosca Structure 2020 (Dassault Systèmes, Vélizy-Villacoublay, France) and the process of minimization of the deviation from a reference stress value was based on the following hypothesis by Neuber: the optimum form of a component is achieved when the stresses running along the considered surface zone is fully constant (stress homogenization). Unlike conventional optimization iterations, our approach eliminated the need for modifying the CAD parametric model in each iteration. Instead, we directly adjusted the coordinates of the design variable nodes, thereby eliminating the time-consuming CAD parametric model modification, and then we repeated finite element pre-processing and the interpolation process for the temperature field.

Furthermore, after optimization by Tosca Structure' shape optimization, we obtained the mesh nodes file of a vent hole, and then the coordinates of the finite element nodes were obtained by using the user-developing languages of the ANSYS 19.2, namely, the APDL scripts. Then, we reconstructed the vent hole boundary in CAD software by combining spline curves with the finite element nodes coordinates from the previous step. In this process, this integration formed the final outcome of our optimization methodology.

**Figure 3.** Flowchart of non-parametric-based shape optimization.

#### **3. Vent Hole of a Gas Turbine Flow Guide Disk**

#### *3.1. Introduction of a Vent Hole*

In this study, we focused on optimizing a gas turbine flow guide disk vent hole, which served as a representative example. Figure 4 illustrates the schematic diagram of a flow guide disk, wherein the vent holes play a crucial role in supplying cooling air to the turbine rotor. It is important to note that the location and configuration of these vent holes significantly impact the performance of the air system, as well as the fatigue life of the disk.

**Figure 4.** Schematic diagram of the gas turbine flow guide disk.

Considering that the positioning and size of the vent holes often result in considerable circumferential and radial stresses, the structural designer devised a solution to mitigate the stress concentration at the hole edges. As depicted in Figure 5, the designer initially employed a design approach featuring a straight line superimposed on multiple circular arcs for the hole, effectively reducing the stress concentration. By addressing these design considerations and optimizing the vent hole configuration, we aimed to enhance the overall performance and durability of the gas turbine system.

**Figure 5.** Schematic diagram of the vent hole (initial design).

#### *3.2. Strength Evaluation of the Initial Design*

The material used for the flow guide disk was the FGH95 Ni-based superalloy, which is commonly used in aerospace and gas turbine applications, especially for aero engines in China. Some of the mechanical properties of the FGH95 superalloy are shown in Table 1.

**Table 1.** Mechanical properties of FGH95.


The flow guide disk had N holes evenly distributed in the circumferential direction, and the structure and loads exhibited cyclic symmetry. To optimize the design process and save on computational time, the calculation model considered a 1/N cyclic symmetric section, which included a vent hole.

Since the vicinity of a hole was the focus of this study, in the finite element model, the mesh near the hole edge was refined based on a mesh sensitivity analysis, as shown in Figure 6. It was observed that the maximum equivalent stress near the hole edge converged when the mesh size was approximately 0.2 mm. Then, the local mesh size for strength evaluation around the hole edge was determined to be 0.2 mm, while a mesh size of 2 mm was used in the remaining regions to reduce the computational analysis time. The finite element analysis model consisted of a total of 182,490 elements, and 289,011 nodes, as depicted in Figure 7. In Figure 7, the x-axis represents the radial direction, the y-axis represents the circumferential direction, and the z-axis represents the axial direction.

The finite element model was subjected to centrifugal and temperature loads. The centrifugal loads were applied to the model in the form of a rotational speed of 20,868 r/min, while the temperature loads were applied to the nodes, and the temperature distribution was as shown in Figure 8. The axial and circumferential displacements of all the nodes on face A, as well as the axial displacements of all the nodes on face B, were constrained. Additionally, the cyclic symmetry constraint was applied to both the lower and upper boundaries.

**Figure 6.** Mesh sensitivity analysis of the finite element model.

**Figure 7.** Finite element mesh of the gas turbine flow guide disk.

**Figure 8.** Temperature distribution of the gas turbine flow guide disk.

Under the given load conditions, the linear elastic finite element simulation analysis was performed with a calculation time of approximately 5 min on an Intel Xeon E5-1620 V4 Core 32 G memory computing platform. The equivalent stress distribution of the flow guide disk was as shown in Figure 9. The maximum equivalent stress was 1513.3 MPa, and it was located at the R3 circular arc section of the hole edge, as shown in Figure 5. This indicated that the edge of the vent hole served as a weak point in the fatigue life of the flow guide disk, and the stress distributions along the upper and lower edges of the hole were closely symmetric with respect to the middle line of the enlarged position shown in Figure 9.

**Figure 9.** Equivalent stress distributions in the gas turbine flow guide disk's initial design.

In order to facilitate the subsequent design optimization, the detailed stress distribution and composition of the hole's front edge were extracted according to the normalized path of 0~0.5~1, as shown in Figure 9. Figure 10 illustrates the stress distribution along the vent hole edge of the initial design, encompassing both the equivalent stress and normal stress. An analysis of the figure revealed the following: the maximum circumferential stress (normal stress on the Y-axis) measured 1392 MPa, while the maximum radial stress (normal stress on the X-axis) amounted to 791 MPa. The circumferential stress distribution closely resembled the equivalent stress and surpassed the radial stress, implying that the circumferential stress served as the predominant stress component. Consequently, to mitigate the maximum stress value at the hole edge, the overall optimization approach proposed in this paper focused on reducing the curvature of the curve at the R3 edge of the hole and employing the "compress the upper and lower edge, stretch the left and right edge" strategy. This study incorporated methods based on geometric parameters and non-parametric techniques.

**Figure 10.** Stress distribution along the vent hole edge in the initial design.

#### **4. Design and Optimization**

#### *4.1. Geometric Parameter-Based Optimization*

After conducting the aforementioned strength analysis, the vent hole scheme was redesigned in this study with the aim of reducing the stress value at the hole edge. As illustrated in Figure 11, the vent hole was modified to an elliptical shape based on our engineering experience, with the semi major axis denoted as '*a*' and the semi minor axis denoted as '*b*'. The major axis of the ellipse defined the circumferential direction of the vent hole, while the minor axis determined the radial direction. Since the air system imposed restrictions on the vent hole cross-sectional area, the cross-sectional area of the vent hole remained constant at '*S*'. This modification addressed the stress-related concerns and ensured compliance with the vent hole cross-sectional area requirements.

**Figure 11.** Schematic diagram of the vent hole based on the geometric parameters for vent holes.

In this study, the mathematical model for geometric parameter-based shape optimization can be represented by Equation (1), as follows:

$$\begin{array}{l}\min \ \sigma\_{\epsilon\eta\upsilon,\text{max}}\\hline\ a\\s.t.\ \mathcal{S}=\mathcal{S}\_{0}\end{array}\tag{1}$$

In Equation (1), *σeqv*,max is the maximum equivalent stress at the hole edge, *a* is the semimajor axis of the ellipse, and *S* represents the cross-sectional area of vent hole, which is equal to a specific value, denoted as '*S*0', according to the requirements of the air system on the vent hole. In this case, *S*<sup>0</sup> was equal to 84.198 mm2. The shape optimization method based on the geometric parameters necessitated adjustments to the relevant geometric model parameters, CAD model re-meshing, and interpolation of the temperature field during each iterative sub-step of the optimization process. Because the cross-sectional area remained basically unchanged, the variations in the vent hole had no significant impacts on the temperature field, and therefore, the temperature load of the finite element model remained unchanged during the optimization iterations. Given that the number of design variables was limited to only one, the NLPQL gradient algorithm [30] was employed in this study to expedite the optimization calculations. Figure 12 illustrates the comparison between the initial design scheme and the optimal design scheme achieved through geometric parameter-based optimization. Furthermore, Figure 13 depicts the distribution of the equivalent stress in the flow guide disk and the vent hole when employing the optimal design scheme based on the geometric parameters.

**Figure 12.** Comparison of the schemes between the initial design and the geometric parameter-based optimal design.

**Figure 13.** Equivalent stress distribution in the geometric parameter-based optimal design.

Table 2 presents a comparison of the optimization results between the initial design and the design achieved through geometric parameter-based optimization. The optimization results indicated that the semi-major axis (*a*) of the ellipse after optimization was 7.3 mm, corresponding to a semi-minor axis value of 3.67 mm. The area of the vent hole was 84.166 mm2, which was 0.04% smaller than the initial value of 84.198 mm2. This slight difference was acceptable within the engineering requirements for the design dimensions of the vent hole, which are typically specified to two decimal places, and it could be accommodated by the air system. The maximum equivalent stress at the hole edge in the optimal design based on the geometric parameters was 1207 MPa, reflecting a significant reduction of 20.2% compared to the initial design, which recorded a maximum equivalent stress of 1513 MPa. This improvement highlighted the effectiveness of the geometric parameter-based optimization approach in mitigating stress levels at the hole edge.

**Table 2.** Comparison of the results between the initial design and the geometric parameter-based optimal design.


#### *4.2. Non-Parametric-Based Optimization*

For the non-parametric shape optimization of the vent hole, all nodes located along the hole's edge were chosen as design variables. To ensure sufficient mesh quality for the finite element analysis, adaptive mesh technology was employed in Tosca Structure. This technique enables the cells surrounding a design's variable nodes to move adaptively within a specified range, thereby maintaining a high-quality mesh around a vent hole. In this study, the optimization objective was to minimize the maximum equivalent stress at the vent hole's edge in the flow guide disk. The specific mathematical model for this optimization process can be expressed by Equation (2), as follows:

min *σeqv*,max *find Nodei*,*x*,*<sup>y</sup>* (*i* = 1... *j*) *s*.*t*. *S* = *S*<sup>0</sup> (2)

In Equation (2), *σeqv*,max is the maximum equivalent stress at the hole edge, *Nodei*,*x*,*<sup>y</sup>* is the location of node *i* of the hole edge, and *S* is the cross-sectional area of the vent hole, which is expected to be kept the same during the optimization iterations. To meet the air system's restriction requirements for the vent hole's cross-sectional area, the constraints were transformed to maintain the optimized structure volume equal to the initial design volume. Due to the manufacturability requirements of the vent hole, it needed to be aligned with the basic hole scheme in the axial direction. To meet the engineering manufacturing requirements, constraints were applied to the movement of these design variables. The nodes located at the same angular positions along the hole edges were selected as a group. Figure 14 illustrates the configuration of the four node groups (red nodes) established in this study, totaling 66 groups, ensuring the consistent direction and magnitude of movement for the nodes within the same groups. If this constraint was not imposed, the direction and magnitude of movement for each node within the same group would have been inconsistent, and the hole would not have been smooth along the axil direction. As a result, the optimization outcome could not have been manufactured.

**Figure 14.** Equivalent stress distribution of the geometric parameter-based optimal design.

The variations in the equivalent (von-Mises) stress during the iteration process is shown in Figure 15. We can see that after 46 iterations in Tosca Structure, the equivalent (von-Mises) stress of vent hole converged to 1130 MPa, and we successfully obtained the optimal results for the vent hole in the flow guide disk.

**Figure 15.** The variations in the equivalent (von-Mises) stress during the iterative process.

We generated an optimized 3D model by using a self-developed APDL script and spline curves. Figure 16 displays a comparison between the initial design scheme and the optimal design achieved through the non-parametric-based optimization approach. It was evident that the optimized scheme based on non-parametric optimization closely resembled the "ellipse-like" scheme derived from the geometric parameter-based optimization. Furthermore, Figure 17 presents the distribution of the equivalent stress in the flow guide disk and vent hole, highlighting the improvements achieved with the non-parametric-based optimal design.

**Figure 16.** Comparison of the schemes between the initial design and the non-parametric-based optimal design.

Table 3 provides a comprehensive comparison of the optimization results between the initial design and the optimal design achieved through the non-parametric shape optimization approach. The maximum equivalent stress at the hole edge in the non-parameter-based optimal design was reduced by 24.5%, measuring 1142 MPa compared to the initial design's stress of 1513 MPa. The non-parametric-based optimal design showcased the successful reduction in stress concentrations, resulting in a more uniform stress distribution throughout the vent hole. These improvements highlighted the effectiveness of the non-parametric approach in achieving stress reductions and enhancing the overall performance of the design. After optimization, the cross-sectional area of the vent hole was approximately 84.304 mm2, which was approximately 0.13% larger than its initial value of 84.198 mm2. This discrepancy was attributed to the error introduced by using spline curves to describe

the optimized vent hole mesh boundary. However, this could be accepted by the air system requirements.

**Figure 17.** Equivalent stress distribution of the non-parametric-based optimal design.

**Table 3.** Non-parametric shape optimization-based method before and after the optimization results.


#### **5. Discussions**

Two optimization methods, namely, geometric parameter-based optimization and non-parametric optimization, were employed to optimize the vent hole in a gas turbine flow guide disk. Figure 18 illustrates the design schemes of the initial design and the optimal designs obtained through the two different optimization methods. The results demonstrated significant changes in the optimal design schemes compared to the initial design. Notably, the optimal designs generated by both methods exhibited similarities, resembling elliptical shapes with slight variations in curvature at certain locations. This indicated that both optimization methods converged towards a similar optimal design solution, emphasizing the effectiveness of both approaches in achieving the desired design modifications.

Tables 2 and 3 present the results obtained using the two different optimization methods utilized in this study. The following observations could be made:

(1) Both optimization methods resulted in an optimal vent hole design that exhibited a more uniform stress distribution compared to the initial design.

(2) Both optimal designs effectively reduced stress concentrations. The geometric parameter-based optimal design achieved a 20.2% reduction in the maximum equivalent stress, while the non-parametric-based optimal design achieved a greater reduction of 24.5% compared to the initial design. It was worth noting that the non-parametric-based shape optimization method proved to be more efficient in optimizing the maximum equivalent stress at the hole edge.

**Figure 18.** Comparison of the schemes between the initial design and two different optimization methods.

These findings demonstrated the effectiveness of both optimization methods in improving stress distribution and reducing stress concentration, with the non-parametric-based method showcasing a higher efficiency in stress optimization.

In order to gain a more detailed understanding of the differences between the two optimization methods, Figure 19 provides the equivalent stress distribution along the front edge of the vent hole for the initial design and the two optimal designs, and the following observations can be made:

(1) In comparison to the elliptical scheme, the non-parametric-based optimization resulted in a smaller curvature at the location where the maximum equivalent stress occurred. This led to a more uniform stress distribution and a lower maximum equivalent stress value.

(2) While the geometric parameter-based design achieved a 20.2% reduction in the maximum equivalent stress compared to the initial design, it is important to note that the establishment of the optimized hole shape boundary relied heavily on the engineering experience of the designer. Different designers may utilize different hole boundaries, which can result in varying design optimization outcomes.

**Figure 19.** Comparison of the equivalent stress distribution along the vent hole edge between the different designs.

These findings highlighted the advantages of the non-parametric-based optimization method as it provided greater flexibility in achieving an optimal stress reduction and uniform stress distribution while mitigating the dependence on individual designer experience.

Based on the fundamental knowledge of fatigue theory, it is known that reducing the maximum equivalent stress and enabling a more balanced stress distribution along the vent hole edge through optimization is beneficial for improving the fatigue life of a vent hole. A comparison of fatigue life of the vent holes before and after optimization is presented in Table 4, and it is based on our institute's fatigue database for the FGH95 Ni-based superalloy. It can be observed that the geometric parameter-based optimal design resulted in an increase in fatigue life from 3701 cycles to 11,066 cycles, representing a 199% improvement, while the non-parametric-based optimal design yielded a fatigue life of 14,486 cycles, indicating a 291% improvement.

**Table 4.** Comparison of fatigue life between the different designs.


#### **6. Conclusions**

This study established design optimization models and processes for vent holes using geometric parameter-based and non-parameter-based designs. Two shape optimization methods have been presented, and the differences between the two optimal designs have been studied. The key conclusions from this study are as follows:

(1) The stress concentration in the vent hole was reduced and the stress distribution became more uniform after optimization compared to the initial design. This optimization method can be used as a reference for the design optimization of similar engineering structures.

(2) The geometric parameter-based optimal design achieved a 20.2% reduction in maximum equivalent stress compared to the initial design, while the non-parameter-based optimal design achieved a 24.5% reduction.

(3) The non-parametric-based shape optimization method demonstrated better optimization results and effectively avoided the poor design optimization results caused by the lack of relevant structural and strength design experience among designers.

(4) Through the application of the two optimization methods, the geometric parameterbased approach resulted in an improvement in the fatigue life of the vent hole from 3701 cycles to 11,066 cycles. Additionally, the non-parametric-based shape optimization yielded a fatigue life of 14,486 cycles.

In conclusion, this study has successfully established effective design optimization models and processes for vent holes and has provided valuable insights into the differences between geometric parameter-based and non-parameter-based optimization methods. These findings will be useful for future engineering designs.

**Author Contributions:** Conceptualization, C.W., Y.Z. and D.M.; methodology, C.W., Y.Z. and H.Z.; investigation, Y.Z. and Z.Q.; writing—original draft preparation, C.W., Y.Z., D.M. and Z.Q.; writing review and editing, C.W., Y.Z., D.M., Z.Q. and H.Z.; supervision, C.W. and D.M.; project administration, D.M.; funding acquisition, D.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Independent Innovation Special Project of AECC (grant number: ZZCX-2018-017).

**Data Availability Statement:** Not applicable.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Luis Fabian Urrego 1, Olimpo García-Beltrán 2, Nelson Arzola <sup>3</sup> and Oscar Araque 1,\***


**Abstract:** Aluminium alloy (AA2024-T4) is a material commonly used in the aerospace industry, where it forms part of the fuselage of aircraft and spacecraft thanks to its good machinability and strength/weight ratio. These characteristics allowed it to be applied in the construction of the structure of a pilot plant to produce biological extracts and nano-encapsulated bioproducts for the phytosanitary control of diseases associated with microorganisms in crops of *Theobroma cacao* L. (Cacao). The mechanical design of the bolted support joints for this structure implies knowing the performance under fatigue conditions of the AA2024-T4 material since the use of bolts entails the placement of circular stress concentrators in the AA2024-T4 sheet. The geometric correction constant (*Y*) is a dimensionless numerical scalar used to correct the stress intensity factor (SIF) at the crack tip during propagation. This factor allows the stress concentration to be modified as a function of the specimen dimensions. In this work, four compact tension specimens were modeled in AA2024-T4, and each one was modified by introducing a second circular stress concentrator varying its size between 15 mm, 20 mm, 25 mm, and 30 mm, respectively. Applying a cyclic load of 1000N, a load ratio R=-1 and a computational model with tetrahedral elements, it was determined that the highest SIF corresponds to the specimen with a 30 mm concentrator with a value close to 460 MPa.mm0.5. Where the crack propagation had a maximum length of 53 mm. Using these simulation data, it was possible to process each one and obtain a mathematical model that calculates the geometric correction constant (*Y*). The calculated data using the new model was shown to have a direct relationship with the behavior obtained from the simulation.

**Keywords:** crack; fatigue; geometric factor; support vector regression; pilot plant

#### **1. Introduction**

*Theobroma cacao* (Cocoa) is considered one of the most important raw materials in international trade; it is a source of foreign exchange in 58 producing countries, highlighting that 89% of this production is found in Ivory Coast, Ghana, Indonesia, Nigeria, Brazil, Ecuador, Malaysia and Cameroon [1]. Cocoa production is sometimes affected by environmental, physical, and chemical factors and inadequate pest and disease control [2]. This crop is mainly infected by disease-causing microorganisms, among which *Moniliophthora roreri* and *Phytophthora* spp. stand out, which are the two main risk factors that directly affect annual cocoa production [3,4]. This is why management alternatives for these diseases, such as the use of plant extracts and essential oils (EO), should be sought, being favorable for environmental sustainability and human health [5]. This implies the development of an infrastructure capable of producing the raw material used in the production of bio-based products.

**Citation:** Urrego, L.F.; García-Beltrán, O.; Arzola, N.; Araque, O. Mechanical Fracture of Aluminium Alloy (AA 2024-T4), Used in the Manufacture of a Bioproducts Plant. *Metals* **2023**, *13*, 1134. https:// doi.org/10.3390/met13061134

Academic Editors: Haitao Cui and Qinan Han

Received: 9 May 2023 Revised: 10 June 2023 Accepted: 13 June 2023 Published: 16 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The study of the phenomenon of fracture has its origin in work proposed by Griffith (1921, 1924); the researcher Irwin (1957) made an important advance by proposing the analysis of fracture toughness as a function of stress, the toughness of a material is obtained from the applied stress and the crack length, but given the different test configurations, these values are known as a function of the failure modes called Stress Intensity Factor (SIF) KI (Opening), KII (in-plane shear) and KIII (out-of-plane shear). Estimation of the SIF in materials with linear elastic behavior (LEFM) is possible by quantifying the nominal stress and crack size. The use of the finite element method (FEM) allows testing for various configurations due to the versatility of the method; using the SMART component integrated into the latest versions of ANSYS finite element software, it is possible to perform simulations to verify the behavior of the geometric factor accompanying the determination of the SIF in the different failure modes.

Researchers have worked on this method, Nairn [6] has proposed the analysis of the geometric factor as a correction factor of the general SIF equation and that it can be expressed as a function of crack length (*a*) and width (*W*), on the other hand, for the author Mecholsky, the geometric factor (*Y*) for semi-elliptical cracks in materials of high hardness and brittleness is used to explain the position and shape of the crack because it is a function of an angle (*θ*) between the surface of the crack front and any peripheral point above it [7], otherwise, authors Taylor, Cornetti and Pugno, cataloged this geometric factor not only in terms of geometry but also in terms of the crack notch [8], on the other hand, when analyzing fatigue crack propagation, the geometrical factor is included in an expression known as the initial value of propagation for short cracks (*a*0), which according to the authors B. Atzori, P. Lazzarin and G. Meneghetti, occurs at the point of intersection between the change in realized stress (Δ*σ*) and the different values of toughness (Δ*K*), where the geometrical factor is calculated using a simulation in ANSYS 5.6 software [9].

The authors Smith and Scattergood have analyzed ceramic materials defining toughness as a sum of two different toughnesses: (*Kbend*), which is defined as a toughness that is a function of the stress intensity factor and the residual toughness (*Kresidual*) that results from the residual stress field due to strain, an equation is obtained to determine the value of the toughness (*Kbend*) by an exponential equation, which involves a shape factor and the crack depth [10]. However, when analyzing materials with a higher degree of ductility, the empirical approaches of authors J.C. Newman and I.S. Raju obtained an equally accepted behavior for the determination of toughness [11].

On the other hand, when testing chromium steels, authors Nix and Lindley determined that the behavior of the shape factor (*Cs*) was also exponential in nature, where the basis again was the ratio ( *a*/*c* ), where (*a*) is the crack depth, and (*c*) is the crack length, where the values of this ratio were previously calculated by subjecting chromium steel specimens On the other hand, when testing chromium steels, the authors Nix and Lindley established that the behavior of the form factor, identifying that for the tests developed, the fracture toughness equation must involve a factor called (*Mf*), which is a crack front correction factor [12].

These correction factors have been the product of a rigorous algebraic analysis, where the basis of each of the equations that describe these factors, part of the simulations in finite element programs, which provide the input values necessary to apply the respective analysis algorithms, as proposed by the authors Clarke, Griebsch and Simpson, who explain how it is possible to glimpse different situations of mechanical nature by means of the Support Vector Regression (SVR), algorithm, this algorithm allows from some input values in the Cartesian plane, to obtain an equation that adjusts to the distribution of given values, considering the dispersion (ξ) of the same one, and the margin (ε) between the support vectors [13]. This algorithm uses the principles of Lagrangian optimization, simultaneously involving Kernel functions, which can have a polynomial, Gaussian or Sigmoidal nature, and as explained by researchers Schölkopf and Smola [14], authors Heydari and Choupani, indicate a correction factor for fracture toughness of a logarithmic nature based on the rate of energy release [15], authors El-Desouky and El-Wazery, indicated fracture toughness

for materials with a high degree of brittleness, using a fifth-degree polynomial (*F*1), whose variable is the relation ( *a*/*W* ), where (*a*) is the crack length and (*W*) is the length of the cross-sectional area [16].

The models for the geometric factor for compact specimens that are currently planned in the literature are based on specimens that have a single stress concentrator.

Aluminum alloys are widely used in industry to reduce the weight of structural components in machine design. Their light weight, easy machinability and excellent fatigue strength make these alloys ideal materials for manufacturing in the modern world. One of the most recent applications is in the high-speed rail industry, where 5083P-O Aluminium alloy is involved in the design of train bodies, leading to increased running speed and improved assembly performance [17].

The goal of the present work is to evaluate the behavior of this parameter when there are two stress concentrators arranged in the specimen geometry since it has been visually observed that the crack rupture direction is affected by the size and location of the concentrator with respect to the propagation plane. In the same way, the study of the geometrical factor implies a study of the stress intensity factor at the crack tip, which means that this last parameter must also be analyzed for the geometrical conditions proposed. The main differentiating feature lies in determining the stress intensity factor by applying the FEM method and then processing this data by means of the SVR algorithm. Applied to non-standardized specimens with variable diameter stress concentrators. Allowing the data obtained by the simulation to be used within the Nadaraya–Watson estimator to find a mathematical model that explains the behavior under variable tensile loading and crack propagation.

The importance of the model to be developed consists in its use as an element for joining bolted joints in structures with diameter variation with the purpose of establishing life projection parameters useful for design processes.

#### **2. Materials and Methods**

#### *2.1. Material and Specimens*

In this article, a modified compact tension specimen (MCTE) was used based on the ASTM-E399 standard, where the value of W was set at 100 mm and its thickness at 10 mm. In the same way, a stop hole of variable diameter was set in; in the same way, a stop hole of variable diameter was fixed in 4 different specimens; therefore, 4 simulations were performed for each modified compact specimen. The dimensions are shown in Figure 1. For the configuration of the SMART method used in ANSYS, it was necessary to characterize the Aluminium alloy (AA 2024-T4). This is hardened by a thermal aging process and presents the mechanical properties contained in [18] and presented in Table 1; similarly, the author Zyad Nawaf Haji [19] characterized the fatigue deformation parameters for AA2024-T4. These parameters are determined from the Ramberg–Osgood and Coffin–Mason equations that describe well the cyclic behavior of the material, but they are not physical laws [20]. These parameters are presented in Table 2.

**Table 1.** Mechanical properties AA2024-T4 [18].


The mechanical and fatigue properties of the AA2024-T4 used for the ANSYS simulation are shown below.

The authors B.M. Faisal, A.T. Abass and A.F. Hammadi [21] subjected to fatigue tests several specimens of AA2024-T4, determining its characteristic S-N curve, which is presented in Figure 2. Table 3 shows the values for the C and m constants for AA2024-T4 determined by authors Yang Guang, Gao Zengliang, Xu Feng and Wang Xiaogui [22] by subjecting eleven specimens to fatigue tests.

**Figure 1.** Dimensions of modified compact specimens used.

**Table 2.** Fatigue Parameters AA2024-T4 [20].


**Figure 2.** Characteristic fatigue behavior of the material AA 2024-T4 [21].

**Table 3.** Paris' Law Constants [22].


All modified compact tension specimens were subjected to failure mode 1 (opening), which causes a unidirectional cyclic stress perpendicular to the crack plane; using a force of 1000 N, using a Stress Ratio (*R*) = −1, since according to author C.M. Hudson fatigue cracks in Aluminium alloy AA 2024 propagate at a faster rate with *R* = −1 than with *R* = 0 when the same load was applied in both tests. Apparently, the compression portion of the loading cycle accelerates crack growth in this material [23]. Therefore, we proceeded to perform the simulation using ANSYS Mechanical ADPL 21R1, from which we calculated the *KI* values and the crack extension values (*a*) in each sub-step of the solution. It was taken from Equation (1), formulated by R.P. Wei [24] for high-strength aluminum subjected to cyclic axial loading.

$$
\Delta K\_{I\ max} = \Upsilon \Delta \sigma \sqrt{\pi a} \tag{1}
$$

Expressing Δ*KI* as (*KI max* − *KI min*) and Δ*σ* as (*σmax* − *σmin*), it is possible to derive Equation (2) algebraically.

$$K\_{I\max} = Y\sigma\_{\max}(1-R)\sqrt{\pi a} \tag{2}$$

Subtracting for *Y* gives:

$$Y = \frac{K\_{I\text{ max}}}{\sigma\_{\text{max}}(1 - R)\sqrt{\pi a}}\tag{3}$$

Which is a dimensional expression that, for the purposes of the investigation, will be compared with the relative crack length, which is expressed by a ratio *a*/*W* and is also dimensionless.

#### *2.2. Computational Model Using SMART Method*

To achieve a better analysis of the SIF at the crack tip, the meshing was defined using a Patch Conforming Method defined by Tetrahedral elements, these elements, according to the ANSYS usage guide [25], are unique to the SMART method. Tetrahedral elements are 3-dimensional in nature, have a quadratic displacement behavior and are well suited for modeling irregular meshes. The element is defined by 10 nodes that have 3 degrees of freedom in each of the *x*, *y*, and *z* nodal directions see Figure 3.

**Figure 3.** Tetrahedral elements.

The element has plasticity, hyperelasticity, creep, tensile stiffness, good deflection, and tension capabilities. It also has mixed formulation capabilities to simulate deformations of nearly incompressible elastoplastic materials and fully incompressible hyperelastic materials [26]. Similarly, the Refinement method was used to perform sectorized refinements in the crack front and its adjacent areas, as well as in the areas surrounding the circular stress concentrators (stop holes). The grids, the number of elements and nodes of each computational model used are shown below.

The Smart Crack Growth component used in the current research does not require the placement of a pre-crack geometry; however, it does require the involvement of 3 parameters defined within the Fracture-Premesh Crack block that are defined based on the geometry of the notch placed in the specimen. These parameters must be nodal surfaces adjacent to the crack front. The first is defined as a function of the nodes located on the edge of the crack front, which in the present work, we call the front, and interprets the place at which the crack propagation would start; the second, which we call the top, interprets the nodes located on the upper adjacent surface to the crack front; the third, which we call the bottom, interprets the nodes located on the lower adjacent part of the crack front.

Similarly, Figure 4a illustrates the nodes on the crack front (front), Figure 4b illustrates the nodes on the top surface adjacent to the crack front (top), Figure 4c illustrates the nodes on the bottom surface adjacent to the crack front (bottom) and Figure 4d illustrates the location of the pre-crack conditions with the coordinate axes, inferring that the propagation will be along the x-axis.

**Figure 4.** (**a**) illustrates the nodes on the crack front (front), (**b**) illustrates the nodes on the top surface adjacent to the crack front (top), (**c**) illustrates the nodes on the bottom surface adjacent to the crack front (bottom) and (**d**) illustrates the location of the pre-crack conditions with the coordinate axes, inferring that the propagation will be along the x-axis.

In Figures 5–8, the working specimen is presented, showing the number of elements and nodes for each of the variations in the size of the diameter of the stress concentrator elements.

**Figure 5.** MTCE-1 mesh with 101742 elements and 151424 nodes.

**Figure 6.** MTCE-2 Mesh with 104899 elements and 156293 nodes.

**Figure 7.** MTCE-3 Mesh with 107059 elements and 159635 nodes.

The numerical model was defined using the patch-forming method with tetrahedral elements. A general refinement was applied to the entire geometry of each specimen by setting the size of each element to 2.5 mm, and then a surface refinement was applied to each of the concentrators. Within the analysis setup, a time-defined solution was established with 10 s as the time limit and 0.5 s as the time in each solution step.

**Figure 8.** MTCE-4 Mesh with 99726 elements and 149780 nodes.

#### *2.3. Data Processing Using Support Vector Regression (SVR) and Nadaraya-Watson Estimator (NWE)*

The SVR algorithm is an algorithm for linear or non-linear regression of points in the Cartesian plane, whose purpose is to find the equation of the hyper-plane that interpolates all the points, based on the use of Kernel functions *K*(*x*) (see Table 4). Furthermore, considering the margin between vectors (*ε*), the dispersion that exists between the margins and the points furthest from it (*ζ*) and the arithmetic mean of the data (*μ*). These are variables that can be imposed during regressions in this way. The quality of the regressions performed can be measured in terms of the coefficient of determination (*R*2), since the closer this value is to 1, the higher the quality of the regression. The rationale of the SVR algorithm lies in the quadratic optimization of a Gram matrix, which is maximized subject to a real domain condition to find the Lagrange multipliers [26].



The calculation of the weight values *wn* that are calculated by means of the Lagrange multipliers were obtained by means of a code in Python language and from which the approximation was carried out by means of the Nadaraya-Watson estimator. To determine which of the Kernel equations shows similar behavior.

According to author Larroca, F., the Nadaraya–Watson estimator is a type of nonparametric estimator that uses an equation, *m*ˆ(*x*), and a fit parameter, *εo*, to interpolate a series of data in the Cartesian plane, this equation, *m*ˆ(*x*), is given as a function of a weighted average of a Kernel density equation (see Table 4), which is assigned to each of the points taken [27]. However, the author Cai. Z. [28] proposes that the function *m*ˆ(*x*) is altered with weight values, *wn*, to improve its fit, therefore, its general form is given by:

$$Y\_{pred} = \mathfrak{m}(\mathfrak{x}) \pm \varepsilon\_o \tag{4}$$

According to the authors Demir, S. and Toktamis, O. [29], the use of the Kernel within ENW implies changing the dot product between vectors by the expression *xn*−*<sup>x</sup> <sup>h</sup>* , then, the function *m*ˆ(*x*) is posed as:

$$\mathfrak{M}(\mathfrak{x}) = \frac{\sum\_{i=1}^{n} \mathfrak{w}\_{n} \cdot \mathrm{K}\left(\frac{\chi\_{\mathfrak{n}} - \chi}{h}\right) y\_{\mathfrak{n}}}{\sum\_{i=1}^{n} \mathfrak{w}\_{n} \cdot \mathrm{K}\left(\frac{\chi\_{\mathfrak{n}} - \chi}{h}\right)} \tag{5}$$

According to the author Fan, J., *K* symbolizes a Kernel density equation, *xn* is the x-axis coordinate of each point in the Cartesian plane, *yn* is the y-axis coordinate of each point in the Cartesian plane, *n* is the number of points in the plane, and *h* is a parameter known as Bandwidth that controls the accuracy of the interpolating equation with respect to the Cartesian plane data; this last parameter is considered optimal when the highest accuracy is obtained, that is, the smallest mean square error (MSE). This error is calculated as [30].

$$MSE = \frac{1}{n} \left( Y\_{pred} - y\_n \right)^2 \tag{6}$$

Then, the equation interpolating the values is expressed as follows:

$$\mathcal{Y}\_{pred} = \frac{\sum\_{i=1}^{n} \mathfrak{w}\_{n} \cdot \mathcal{K}\left(\frac{\mathbf{x}\_{n} - \mathbf{x}}{h}\right) y\_{n}}{\sum\_{i=1}^{n} \mathfrak{w}\_{n} \cdot \mathcal{K}\left(\frac{\mathbf{x}\_{n} - \mathbf{x}}{h}\right)} \pm \varepsilon\_{o} \tag{7}$$

The analysis procedure used to process the data obtained by using the finite element method is illustrated in Figure 9, where it is shown that the desired Kernel equation and the data obtained from the simulation are used to complete the Gram matrix. This same data will then be used in the NWE to multiply the vector weights obtained from the Lagrangian optimization problem. From this multiplication, a mathematical model is finally obtained. It should be noted that the definition of the dependent and independent variables is maintained from the beginning of the processing and does not change at the end of the processing.

**Figure 9.** Methodology of analysis applied to FEM data.

#### **3. Results**

Figure 10 illustrates the behavior of the SIF in failure mode I with respect to crack length (a):

**Figure 10.** Mode I SIF with respect to crack length (a) for each specimen tested.

Based on the results illustrated in Figure 11, it is observed that MTCE-1 does not present a failure in the concentrator as observed in the other specimens. Therefore, the propagation length was higher in this specimen; this can be explained by the function of a greater cross-sectional area existing in MTCE-1 since the diameter of its circular concentrator is 15 mm, which is the smallest among the other specimens; Figure 11 shows the crack propagation obtained in each specimen.

**Figure 11.** Crack propagation obtained in each specimen, (**I**) illustrates the solution in the first sub-step of the solution, (**II**) illustrates the propagation in the middle sub-step of the solution and (**III**) illustrates the propagation in the last sub-step of the solution.

Figure 12 shows the behavior of the geometric correction factor (*Y*) as a function of the relative crack length (*a*/*W*).

**Figure 12.** Geometric factor compared to the relative crack length.

From these calculated data, a regression was developed by applying the concepts of the SVR and the Nadaraya-Watson estimator (NWE), applying the Epanechnikov Kernel since, according to the authors Chu, C. Y., Henderson, D. J., and Parmeter, C. F, this is the equation that shows the highest efficiency when interpolating data placed in a Cartesian plane [30].

This model was found using a value for Bandwidth (*h*) equal to 0.47 and a value for *ε<sup>o</sup>* of 0.4, whose associated curve is illustrated in Figure 13.

**Figure 13.** Mathematical model relating the geometric correction factor as a function of relative crack length.

This model has an MSE of 11.2%, which may be due to the high dispersion presented when 0.2 ≤ *a*/*W* ≤ 0.4 as from this value, in Figure 12. A high dispersion of the data is observed; this dispersion is caused by the variation in the diameter of the circular stress concentrator used in each MTCE analyzed, so it is prudent to infer that when the specimens analyzed present unique geometric characteristics, the mathematical models obtained must be similarly unique for each MTCE, since for the specimens that are standardized, there are already validated models for their geometries.

#### **4. Discussion**

With the model found in Figure 13, the design of a bolted joint can be proposed from Equation (2), where the hole through which the bolt will pass is assumed to be a circular stress concentrator whose diameter is among those used in each of the MTCE specimens. Then. Starting from Equation (2), we have that:

$$K\_{I\max} = Y\sigma\_{\max}(1-R)\sqrt{\pi a} \tag{8}$$

By subtracting for *σmax*, the expression is left as:

$$
\sigma\_{\text{max}} = \frac{K\_{I \text{ max}}}{Y(1 - R)\sqrt{\pi a}} \tag{9}
$$

Figure 14a illustrates one of the 10 mm diameter holes used in the assembly of the pilot plant structure, and Figure 14b illustrates the force direction (black) in the hole due to the estimated weight of the equipment. The load value used is 1000 N, corresponding to the force exerted on the joint under cyclic loading with R = −1.

**Figure 14.** (**a**) support bolt hole, (**b**) forces assumed for joint design.

Assuming that the material of the structure will be AA-2024-T4, then the joint dimensions can be calculated by the following expression:

$$A = \frac{\left(\Upsilon\_{Al}\right)\left(\frac{-6\left(\frac{a}{w}\right)^2 + 2.75\left(\frac{a}{W}\right) + 0.56}{-3.39\left(\frac{a}{W}\right)^2 + 0.012\left(\frac{a}{W}\right) + 1} + 0.4\right)(1 - R)\sqrt{\pi a}}\tag{10}$$

By performing the calculations using Equation (10), we obtain that the fatigue resistance area of the joint should not be less than 76.2 MPa. Figure 15 illustrates a rendering of the structure model designed.

It is observed that both the behavior of the SIF and the crack length (a) found in this research have a similarity with the data found by the authors Alshoaibi, Abdulnaser M. [31], in this study, the behavior of crack propagation is analyzed in aluminum alloy 7075-T6 specimens, which were subjected to a force of 20 KN and a stress ratio of 0.1, finding crack lengths of almost 30 mm, and with values for the SIF of about 7000 KPa.mm0.5.

The reason for the difference between this study and the current one lies firstly in the composition of each of the alloys, the 7075-T6 alloy being stronger, and secondly, the high magnitude of the applied force with respect to the force proposed in the present investigation, which was 1KN, however, the behavior of the SIF is like that found in this article. Figure 16 illustrates a comparison of the graphs of SIF vs. crack length (a) between the works of the authors Alshoaibi, Abdulnaser M. (a) and the present one (b).

**Figure 15.** Model of the designed structure.

**Figure 16.** SIF data comparison. (**a**) Illustrates a comparison of the graphs of SIF vs. crack length between the works of the authors Alshoaibi, Abdulnaser M. and (**b**) the present one.

Looking at the scale of the values obtained in Figure 16. It is possible to compare them with the results found by the authors J. M. D. Rahmatabadi, M. Pahlavani, A. Bayati and R. Hashemi [32], who, in their work, propose the use of a standardized compact stress specimen with a size between 23.75 mm × 22.8 mm using a dual phase melt at 770 ◦C of Mg LZ71 and Mg LZ91 alloys subjected to quasi-static loading which varied between 0 N and 900N. In the results, the authors show that for a load of around 900 N a *KIC* = 18.9 MPa·m1/2 is generated. If we look at Figure 16, we have that for all the studied MTCEs, the approximate *KIC* value is 70 MPa·mm1/2, which dimensionally equals 2.21 MPa·m1/2. This clear difference between the *KIC* values could be explained by the thickness of the plates used, being the thinner plate used by the authors J. M. D. Rahmatabadi, M. Pahlavani, A. Bayati and R. Hashemi, which is 1.7 mm thicker than the one with the smallest stress cross-sectional area, while the one used in this work is 10 mm thick, which implies that much less stress will be concentrated on the 10 mm plate.

#### **5. Conclusions**

Computational models were generated for each proposed MTCE using a high density of elements in the meshes, with an average of 103356 elements. Likewise, the use of sectored refinements in the notch and in the circular stress concentrator facilitated the analysis of each of the models. Observing the behavior of the SIF in each MTCE, it can be deduced analytically that the presence of a circular stress concentrator alters the area where the stress applied to the specimen affects and therefore alters the SIF. This coincides with the behavior observed in Figure 9 since the presence of a smaller diameter concentrator implies a larger area of incidence of the cyclic stress; therefore, if there is a larger area of incidence, the stress along the cross-sectional area is smaller, and therefore, the SIF takes longer to reach high values, which results in a longer propagation length. When determining the geometric factor *Y*, it was observed that it increases as the diameter of the circular stress concentrator also increases; this is closely related to the behavior of the SIF since this factor was determined based on the data obtained from the SIF in failure mode I. The mathematical model found presented an MSE of 11.2%, fitting the geometric factor (*Y*) with a function that depends on the relative crack length (*a*/*W*).

From the propagation instants in Figure 10, it can be observed that the direction in which the crack propagates changes as the concentrator hole becomes larger. Based on this idea, it can be deduced that as the radius of the hole is farther away from the crack front plane, the propagation angle tends to maintain values close to 0◦, while if the distance between the hole radius and the midline of the specimen becomes shorter, the propagation angles tend to be much higher.

It would be important to analyze the behavior of the propagation angle direction under the same geometrical conditions of the specimens shown in the present work, as well as to establish a relationship between the size of the second concentrator and the crack propagation direction.

The behavior of the crack propagation direction is determined by the crack length. If it is observed in Figure 10 (III) for each of the specimens, the propagation angle with respect to a horizontal axis becomes tighter as the crack length tends to be smaller. This behavior could be explained as a function of the reduction of the cross-sectional area because of the diameter change in the second circular concentrator.

The joint designed using Equation (10) for the bioplant structure proved to be useful at the time of construction since the structure will be subjected to vibration, which produces alternating stresses that fall on the bolted joint. The nomenclature used in the manuscript is shown in Table 5.


**Table 5.** Nomenclature.

**Author Contributions:** Conceptualization, L.F.U., O.G.-B., N.A. and O.A.; Formal analysis, L.F.U., O.G.-B., N.A. and O.A.; Investigation, O.G.-B., N.A. and O.A.; Methodology, L.F.U., O.G.-B., N.A., and O.A.; Resources, O.G.-B. and O.A.; Writing—original draft, L.F.U.; Writing—review&editing, O.G.-B., N.A. and O.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the funding provided for this research from the Department of Research of the University of Ibagué, the Ministry of Science, Technology and Innovation, the Ministry of Education, the Ministry of Industry, Commerce and Tourism, and ICE-TEX, Program Ecosistema Científico-Colombia Científica, from the Francisco José de Caldas Fund, Grand RCFP44842-212-20184.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank the University of Ibague research department, the funding provided for this research from the Ministry of Science, Technology and Innovation, the Ministry of Education, the Ministry of Industry, Commerce and Tourism, and ICE-TEX, Program Ecosistema Científico-Colombia Científica, from the Francisco José de Caldas Fund, Grand RCFP44842-212-20184.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
