## 1 Introduction
### Problem Definition  
A multilinear material model acts as a piece-wise linear approximation of the continuous material curve. Using a more densely selected set of points ensures better alignment to the original curve and, therefore, better accuracy of the resulting stress values obtained from FEA. Especially when modeling components with complex material properties, such as surface-hardened steel components, balancing accuracy with computational efficiency can become important since more complex models take a longer time to complete. Figure 1: An illustrative example of a multilinear curve used for defining a multilinear material model. Surface hardening is commonly employed in highly stressed parts of components to enhance their durability and performance. Critical regions are selectively heat-treated to increase their hardness while preserving the integrity and ductility of the rest of the components [11]. For example, metal gear teeth often undergo surface hardening, as they experience cyclic loading that can lead to material fatigue and gear failure [12]. Such materials, possessing gradually varying material properties, are often referred to as functionally graded materials (FGMs) [13]. Surface-hardened gears and surface-hardened components in general exhibit varying levels of hardness, strength, and residual stresses from the surface towards the core, with the exact distribution requiring precise measurement and experimental testing. However, due to the complexities and high associated costs of conducting cyclic experiments, the common approach is to estimate cyclic and fatigue properties from tensile test data or hardness measurements [14; 15; 16]. Furthermore, the utilization of empirical expressions [17; 18] can provide estimations of the hardness distribution across the surface-hardened component, from its outer layer to its core, based on the heat treatment parameters applied during the surface-hardening process. Employing the mentioned estimation techniques to establish correlations between hardness and other material properties enables the comprehensive characterization of the material of such components, which can then be effectively modeled to predict their mechanical behavior. Significant advances in determining the mechanical behavior of surface-hardened components have been made, from early analytical methods [13; 19] to more sophisticated computer-aided numerical methods [12; 20]. A widely adopted strategy for the numerical modeling of surface-hardened components, or, more generally, functionally graded materials, is the multilayer method [12; 20; 21; 22], which involves dividing geometry into layers, with each layer characterized by constant homogenous properties. Figure 2 illustrates a typical distribution of Brinell hardness observed in surface-hardened gears with the discretized hardness profile that is then used in FEA to apply the multilayer method. Regarding the multilayer method, [PERSON] and [PERSON] [20] introduced a finite element model to predict the monotonic and cyclic behavior of steel specimens under axial loading. Two models are proposed within the study: a two-layer model and a four-layer model. In the two-layer model, the surface-hardened specimen is divided into a hardened surface and a core, with each layer modeled as a homogeneous material having its own stress-strain curve derived from experimental tests. The four-layer model incorporates two additional layers featuring homogeneous material properties, with each layer characterized by its own stress-strain curves obtained through the interpolation of stress-strain data from the two-layer model. Similar to the two-layer model mentioned earlier, [PERSON] [21] proposed a model to consider the effects of surface hardening in the fatigue strength assessment. The surface layer of the component is modeled as homogenous and exhibits high hardness with an accompanying compressive residual stress state.