Comparative Analysis of Fatigue Energy Characteristics of S355J2 Steel Subjected to Multi-Axis Loads

In this study, a novel test system for estimating bending and torsion fatigue under selectable kinematic and dynamic loading modes was constructed. Using S355J2 steel specimens, a series of tests were conducted to determine material sensitivity to different load paths and loading modes. The experimental results were supplemented with the results of numerical analyses, on the basis of which the components of strain and stress tensors for subsequent analyses were determined in the entire working part of the specimen. An original method for determining specific strain energy components was used. The experimental results showed the grouping of data according to the mode of loading chosen. This could signify that the selected fatigue models are sensitive to certain loading scenarios. We performed in-depth data analysis and complex numerical simulations, formulating likely explanations for the observed effect.


Introduction
Non-monotonic alternation of stress states may lead to permanent changes in the structure of materials and is often the cause of the limited functionality of components. The vast majority of machines and engineering structures are subjected to complex and time-varying operational loads [1][2][3][4][5][6]. Generally, the loads can be dynamic, i.e., caused by forces and moments varying in time. They can also be kinematic loads, caused by time-varying displacements and rotations. In both cases, a complex field of strain and stress arises in the structure. Industrial examples, where such loading scenarios are important, include heavy vehicles used for landscaping, such as diggers or crushers. Operation of such equipment requires completing target kinematic displacement against a uncontrollable and unpredictable dynamic reaction from soil or processed material. A reverse case example can be found in many torsionally loaded shafts, where the kinematic response is stochastic, such as large bore drills in the oil and gas industry. Under the influence of non-proportional loadings, the principal stress and strain directions change cyclically [7,8] and, as proven by various studies [7,9], the fatigue life of a given structure can be significantly reduced. A number of studies on a certain group of materials [10] showed a clear impact of material non-proportional hardening sensitivity on the observed fatigue life. Most experimental studies conducted under multiaxial fatigue load conditions were performed for independent tension: compression and torsion [11,12]. The results of material testing under other loading conditions, for example, for bending with torsion, appear much less frequently in the literature. In most cases, the tests were performed phase-aligned with controlled values of bending and torsional moments amplitudes (dynamic load) [13][14][15][16][17]. Testing of materials subjected to bending and torsional loads, with controlled kinematic load, is much less common; however, some studies [18,19] analyzed essentially aligned kinematic loads causing bending and torsion. Literature sources on material testing in both possible load control variants (dynamic and kinematic) is lacking. In this context, the knowledge of fatigue of materials subjected to phase-dependent kinematic or dynamic loading seems to be underdeveloped, especially given potential significance to many structures and mechanisms. As a part of a research program, material specimens were tested with controlled, with independent kinematic or dynamic loads. A custom-built fatigue testing machine was implemented to facilitate this project. Similar to typical tensile-compression tests with torsion, fatigue tests were performed for both out-of-phase loads and in-phase loads for bending and torsion alone and combinations thereof.
The main objective of this study was to examine the fatigue properties of S355J2 steel using an energetic description. In this work, fatigue tests were conducted on specimens in kinematic and dynamic scenarios with four paths of loadings. The cyclic stress-strain response of the fatigue-tested specimens were analyzed using finite element method (FEM) analysis with the Chaboche model of cyclic plasticity. We determined the effectiveness of the most popular energy models for many variants of cyclic loads based on the original method of determining specific strain energy components to describe fatigue life. The fatigue lives of these specimens were compared between kinematic and dynamic loading test conditions.

Test Specimen
Steel grade S355J2 was supplied for testing in the form of drawn bars with a diameter φ = 16 mm. The chemical composition of this material was previously published [16], whereas its relevant mechanical properties are listed in Table 1. The test specimen is shown in Figure 1. The chosen shape of specimens ensured that the maximum normal and tangential stresses occurred in the smallest cross-section area, making it the most probable location of fatigue crack initiation.

Experimental Procedure
The custom-built test stand we used in oue research is presented in Figure 2. The fatigue machine rests on a desktop base plate, supported by a frame structure. The test specimen is gripped between a tool holder and a column clamp. Each loading lever (bending and torsional) is connected to the holder tool and through a swivel joint, linked with an actuator that provides excitation.  A schematic diagram of the control and conditioning system used in the experiments is shown in Figure 3. It is based in NI cDAQ hardware (National Instruments, Austin, TX, USA).
A test specimen ( Figure 1) was mounted in a tool holder that allows bending and torsional loads to be applied independent of each other. A computerized control system enabled the implementation of any load path and two work modes: kinematic and dynamic. In the former, the loads are delivered as displacements, causing independently controllable torsion X T or bending X B , as well as a any combination of both. The maximum permissible range of lever movements was set to ±9 mm. However, experiments were conducted in the range of ±3 mm. Control and acquisition of displacement values were achieved with dedicated laser displacement sensors, monitoring each load lever. In the second dynamic mode of operation, the specimen was loaded dynamically using independently controlled moment of force: torsion M T and bending M B . In this case, the permissible range of lever displacement was also limited to ±9 mm. The values of bending and torsional moments specified in the experiment were within the range of ±28 Nm. The control and acquisition of the moment values were carried out using appropriately calibrated strain gauges located on the load levers. The specimen neck (its smallest diameter, see Figure 1) caused a concentration of the maximum elastic and plastic strains. The average components of the ε i,j strain tensor at a selected measuring point was measured using resistance strain gauges. However, this required sticking a flat strain gauge on the saddle surface of the specimens' neck. Even with relatively small strain gauges, their deformation and averaged measured deformation values must be considered. We decided to calculate the values of strain tensor components by indirect means, as described in detail below. The control and measurement system also acquired strain value data at the specimen neck. These data were used for preliminary calibration to determine the relationship between lever displacements and strain measured at the neck, which aided in later numerical modelling work. In the conducted experiments, we assumed that the respective ratios of bending and torsion displacements or moments of force amounted to unity. All loading paths used in experiments were fully reversed, as shown in Figure 4, as follows: path I, bending; path II, torsion; path III, proportional in-phase bending with torsion; and path IV, out-of-phase bending with torsion. As the experimental tests were conducted using kinematic and dynamic modes, two different criteria for specimen failure were adopted. For kinematic loads, the number of cycles to specimen failure was defined as the durability corresponding to a 20% change in the measured bending or torsional moment. For dynamic loads, the fatigue life was defined as number of cycles until a 20% change in the displacement of the bending or torsional lever was achieved. All tests were conducted at an excitation frequency of 20 Hz.   Table 2 summarizes the obtained values of fatigue life for specimens tested under controlled dynamic or kinematic loads. Specimens marked as H were subjected to kinematic loading, whereas specimens marked as P were subjected to dynamic loading. For construction materials in an elastic state, stress can be easily calculated directly on the basis of strains (or vice versa, strains on the basis of stress). In this case, the generalized Hooke law can be used. However, this approach ceases to be adequate as soon as the first plastic strain appears. In this case, the relationship between the strain and stress is strictly dependent on the load trajectory. Then, a constitutive model of cyclic plasticity must be used. Constitutive models, despite their diversity in terms of scope of applications and phenomena they embrace, always consist of three basic components.

•
The plasticity condition in which the von Mises equation is the most commonly used where s is the stress deviator, α is the back stress tensor, σ y is the radius of the plasticity surface (the limit of cyclic plasticity in uniaxial tensile test), the mark ":" represents the scalar product of two tensors, and denotes a tensor. • Associated flow rule, which allows the determination of the increase in plastic strain, which, according to Drucker's postulate [20], is normal to the plasticity surface where ds is the stress deviator increase, n is the normal vector for plasticity surfaces n = 3 2 (s−α) σ y , and H is the plasticity module. • The hardening rule. In most cases, to calculate fatigue, the kinematic hardening rule is used, allowing the position of plasticity surface α to be determined after each increment in the stress deviator s. For example, for the non-linear Chaboche model [20,21], the hardening rule takes the following form: where C i and D i (Table 1) are the parameters of the cyclic strain curve function [22] To determine the components of the strain and stress state in specimens prior to experiments, a finite element method simulation was conducted on a numerical model of the test stand, produced in ANSYS 2019 R1 software (Ansys, Canonsburg, PA, USA). This allowed the computation of the field of stress σ(t, x, y, z) and strain ε(t, x, y, z) tensors in the specimen when subjected to loadings in form of moments or displacements.
The numerical model of the fatigue machine consisted of two levers, a holder tool, and a test specimen ( Figure 5). These virtual components were supported and loaded to reflect the kinematics of fatigue machine operation. Although some minor geometrical features of the design were simplified, clamping and loading remained unaffected compared to the actual test stand. A detailed view in Figure 5 depicts the FEM model mesh structure (Table 3) for the virtual specimen, which consisted of 245,561 finite elements distributed over the neck region of the specimen. The finite element mesh was optimized using standard ANSYS procedures. The FEM model uses a diverse finite element mesh that was optimized for the given calculation model.  The cyclic properties and other necessary material data required to introduce the Chaboche constitutive model were determined on the basis of the cyclic strain curves of the S355J2 steel specimens. The loading conditions of the specimen, implemented in experiments under kinematic or dynamic loading, ensured that only normal strains ε yy and shear strains ε xy (and their corresponding stress values) reached significant values (Figures 6 and 7). Other strain and stress values were significantly smaller.  As the chosen shape of the specimen (Figure 1) guided the fatigue fracture location (within the neck region), further analysis only considered the states of the strain and stress in the selected, most stressed point on the surface of the specimen, as indicated by numerical simulations.

Energetic Description of the Experimental Results
The comparison between the fatigue life of specimens tested with two different load modes (kinematic and dynamic) required the selection of the correct fatigue parameter. Among many models found in the literature [22][23][24][25][26], five common energetic models were used: The Smith-Watson-Topper (SWT) ∆W SWT parameter [27] is based on a product of maximum stress σ n,max and principal strain range ∆ε 1 located on the principal strain range plane Chu [28] proposed a fatigue parameter using tensile work σ n,max ∆ε 2 complemented by torsion work τ n,max ∆γ 2 : •

Liu model
For multi-axial loads, Liu [29] proposed a model adopting two different basic modes of fatigue damage. The first is caused by normal stresses, characterized by virtual energy ∆W (I) . The second mode is produced by shear stress, characterized by virtual energy ∆W (I I) . Fracture is expected on a plane where the value of the virtual energy reaches its maximum. The calculation of the ∆W (I) value is preceded by a search for the plane on which the work on the normal strain is the highest. Subsequently, work on shear strain is calculated.
Virtual work is calculated similarly ∆W (I I) . In this case, the value of the work on the shear strain is considered first: Glinka et al. [30] proposed a fatigue parameter calculated as the sum of normal and tangential stress work on corresponding strains, determined by the plane of maximum shear strains: The authors, aiming to include the impact of mean stresses in the parameter in Equation (9), proposed modifications in [31]: where σ f , τ f are the fatigue strength coefficients for tensile compression and torsion, respectively.

Ellyin model
The Ellyin model [32,33] proposes both plastic strain energy ∆W p and positive elastic strain energy ∆W + : For selected cases of proportional and non-proportional loads, the energy is calculated using the following formula: where σ ij and ε p ij are a stress tensor and plastic strain tensor, respectively; σ i and ε e i are the stress and elastic parts of the principal strains, respectively; and H(x) is the Heaviside function.

Calculation and Identification of Specific Strain Energy Components
When applying the Ellyin model, which refers directly to the elasto-plastic strain energy density, it was necessary to calculate it properly using the strain and stress tensor time curves determined previously through numerical simulations.
Previous studies [9,34] presented an application case for the methodology described in [35]. Through the proposed instantaneous power, defined as: a consecutive time sequence of instantaneous steady states can be presented with the following equation: where δW is an increment of work of the internal forces on the infinitesimal elongation/shortening increments. Figure 8 presents the method of calculating the strain energy density for cyclic elasto-plastic loadings. The computational methodology is illustrated on an example of uniaxial tension-compression; however, it is analogous to bending and torsion with respect to the change in the selected strain or stress tensor values. The specified extent of loading variations consists of two ranges: Al-C, which is the primary range, and C-E, which is the complementary range. Within the first, another distinction can be made between the tension relief (A-B) and compression (B-C) ranges. Similar to the primary range, the complementary range also comprises both a compression relief range (C-D) and a tension range (D-E). The area cut off by the power time course, as a result of the work of internal forces during a relief phase, is represented by W I . This is followed by the W I I corresponding to the strain energy density related to the work of external forces throughout the compression phase and the W I I I representing the work of internal forces during the relief after compression phase. Finally, W IV denotes the strain energy density related to the work of the external forces while the material was undergoing the tension phase. Notably, the sum of the areas that ware cut off by of the power time course is not equal for regions located above and below the time axis. By adding W I I and W IV and subtracting the result by the sum of W I and W I I I , the plastic strain energy density can be calculated for the specified pair of ranges (hysteresis loop). This method of identification and calculation of the total strain energy density components [35], as described above, can be applied to the results of the FEM simulations.  Figure 9 depicts an example of a stabilized hysteresis loop and instantaneous power time courses for bending p σ and torsion p τ in the numerical simulations. These were conducted in-phase in kinematic loading mode. Figure 10 shows the corresponding numerical data results from out-of-phase kinematic load simulations. Figure 11 illustrates an example of the instantaneous power time course for bending p σ and torsion p τ , captured in simulations, under dynamic in-phase loading. In addition, Figure 12 presents the corresponding results for dynamic out-of-phase loading simulations.

Discussion
Tables 4 and 5 summarize the results of the numerical calculations of the elastic-plastic strain energy (11) and the fatigue energy parameters from Equations (5) to (10). ASTM standard recommendations were used to assess suitability of the selected fatigue parameter for describing the performed fatigue tests. Following ASTM standard guidelines [36], a linear regression model was adopted (presented in a double logarithmic system) to assess the suitability of the selected fatigue parameter where N f is durability expressed in cycles, W is the value of the strain energy or fatigue energy parameter, and A and B are parameters calculated for this regression model. Table 6 summarizes the calculated values of the model parameters in Equation (16) and the values of the correlation coefficient xy . Values of xy differ depending on the mode used for specimens loading. For kinematic loads, the values are similar and not less than 0.9, with the exception of the Glinka model, where xy = 0.62. For dynamic loads, the values of xy are significantly lower, in the range of 0.47 to 0.8.
All considered energy models achieved good correlation with experimental data regardless of loading path, but only within each loading mode used in the experiments (either kinematic or dynamic). Table 6. The correlation coefficient xy and the regression curve of Equation (16) 9)) and strain energy (Equation (11)) as a function of fatigue life. For all energy models, we noticed that the values obtained in kinematic loading mode were clearly separated from the results acquired in dynamic loading mode. Both groups formed distinct scatter bands in charts, which seemed to be broadly parallel to each other.
All models described above use the work of stress on strains. However, only the Ellyin model considers the hysteresis loop directly in the calculation of strain energy. The SWT, Liu, Chu, and Glinka models propose the use of certain values describing the state of stress and strain that are relatively easy to calculate. However, according to Liu [29], their relationship with the plastic strain energy is abstract. It is especially obvious for non-proportional loads, for which these models were not originally intended for use.
Evaluating the results presented in Tables 5 and 6, we observed that the values of Liu's energy parameters were always greater than strain energy calculated from Ellyin's models. This was noticeable for both loading modes: kinematic and dynamic. For the former, the SWT, Chu, and Glinka energy parameter values were smaller than strain energy calculated from Ellyin's equations. For the latter mode, a different situation occurred. In this case, both SWT and Chu's parameters return values close to Ellyin's, whereas results calculated according to Glinka are noticeably larger.
Comparing the durability under different load methods (kinematic and dynamic), we found that the value of elastic-plastic strain energy or any energy parameter can be significantly different for the same achieved fatigue life.
Searching for the underlying reasons for the difference in strain energy between kinematic and dynamic loading modes, we directly analyzed the time courses of the moment of force and displacements of bending and torsion levers. Figures 16 and 17 depict the total work imputed into the excitation system: levers, holder, and specimen. The rigid design of the test stand ensured that only an insignificant amount of that energy was spent outside of the specimen, dissipated through inter-component clearances, fastened connections, bearings, etc.     3 Figure 17b illustrates that the tests conducted under kinematic load (lines in red and green) returned a near constant value of amplitudes of moments of force for approximately 5/6 of the the fatigue test duration, only to start dropping near the end of the test. The amplitude of displacements ( Figure 17a) did not change throughout the experiment in this mode of loading (kinematic). Throughout the first 5/6 of the test, instantaneous strain energy was constant; however, as the moment of force decreased in the final stage of the experiment, so did the strain energy.
A similar effect was observed for tests under dynamic loading. However, in this case, the observed effects were different. Figure 17a depicts an increase in the amplitude of both lever displacements (black and blue curves) at a constant value of moment loads (Figure 17b). Notably, in both figures, the black lines correspond to the bending and blue to the torsion. Here, at the conclusion of the tests, the amplitude of displacements increased. Since the values of the moment of force amplitude did not change (Figure 17b), the value of the strain energy must increase when dynamic load mode is used.
Among the considered energy fatigue parameters, none consider the temporal transiency of strain tensor and the corresponding stress tensor. As a consequence, the values of strain energy (stress work on strains) were calculated with these models for a given time instance. Thus, the results ignore the changes in the values of the strain and stress tensor components during the fatigue process, affecting the strain energy or any energy parameter calculations.
The numerical constitutive cyclic-plasticity model (Chaboche) used in calculations does not consider an important transiency in material properties caused by damage accumulation or crack initiation, propagation, etc. As a result, the obtained strain and stress tensor values are likely impacted.
A second limitation of the considered energy fatigue models is related to their spatial sensitivity (local effect). This is understood as the need to calculate strain energy only at the most stressed point in the smallest cross-section of the specimen. This, in addition to sampling at an arbitrarily determined time instance, cannot provide complete information on the distribution of strain energy in a representative volume of the specimen. This volume can be considered to be subjected to plastic strains. Figure 18a depict the distribution of mesh nodes at the specimen neck during in-phase kinematic loading. The numerical model expressed the highest levels of stress at the two opposing node regions located on both sides of the neck (marked in red). The plastic stain distribution for this loading case is illustrated in Figure 18b. 10 Figure 17. (a) Change in the fatigue testing machine lever displacement throughout the entire duration of the specimen test. Red and green correspond to the specimen subjected to kinematic loading, whereas black and blue correspond to the fatigue testing machine lever displacement for the specimen subjected to dynamic loading, (b) The moments load the specimen throughout the entire test period. Red and green correspond to the specimen subjected to kinematic loading, whereas black and blue correspond to the specimen subjected to dynamic loading.
In comparison, Figure 19a demonstrates the corresponding distribution of mesh nodes during out-of-phase kinematic loading. Here, the highest stress concentration at the nodes form a ring encircling the narrowest cross-section of the neck. This region corresponds directly to the plastic strain distribution region, as presented in Figure 19b. Similar maximum stress patterns and plastic strain regions were observed for dynamically loaded numerical simulations.
The numerically obtained values of strain energy and energy parameters were expected to be correctly estimated by a single regression model, regardless of loading path and irrespective of loading mode (dynamic or kinematic). Only the former of these assumptions was fulfilled, suggesting that the evaluated common fatigue models might be not entirely applicable in every load case.

Conclusions
• We conducted experimental research of the bending and torsional properties of S355J2 steel under two different modes of load application, which showed the significant limitations of some popular energy parameters used to determine fatigue durability.

•
Despite similarities observed in the fatigue life of specimens tested under kinematic and dynamic modes of load application, neither energy parameters nor strain energy calculations allowed for the common regression model to be implemented for all results. Kinematic and dynamic data required separate regression models.

•
The energy parameters published in the literature were in good agreement with data obtained from different loading paths, but only within specific groups of results obtained from either the kinematic or dynamic mode of loading.
• Among the possible reasons for such behavior, we identified the three most likely causes: 1. The constitutive model of cyclic plasticity (Chaboche) influences the numerically obtained values of strain and stress tensors. The used model assumes that these values are stable and do not change in time, in contrast to actual fatigue damage. This can significantly influence the energy parameter and strain energy values. 2. The described popular fatigue models consider only a singular arbitrary chosen time instance in fatigue life to compute either the energy parameter or strain energy. This results omit the time dependency in strain energy. 3. The spatial sensitivity (local effect) of fatigue parameters proposed by many authors is a limitation affect the accuracy of the results. This is because when calculating many energy parameters as a single parameter, only the most stressed location is considered. As a consequence, certain parts of plastic strain energy deposited in the specimen are omitted.
• The custom-designed and built fatigue test stand will enable further investigations into bending and torsion fatigue behavior of materials under many load paths, including random loads and both kinematic and dynamic load application mode.
As energy models do not consider the change in energy over time associated with the initiation and development of fatigue damage, the variability of energy parameters must be considered during the fatigue test and an algorithm must be developed for estimating energy values over a specific section, area, or volume depending on the analyzed load variant.