1 Introduction

The structural integrity of engineering structures under cyclic loading is an important aspect for different industries. Under cyclic loading conditions, materials experience repeated tensile and compressive load showing different behaviors under monotonic loading [1, 2]. Many researchers examined the behavior of materials and structures which underwent the cyclic tension–compression test. Kim et al. [3] tested the influence of loading rate on the fracture behavior of nuclear piping materials under cyclic loading conditions. In another work, Kim et al. [4] examined the behavior of CF8A austenite cast stainless steel in a tension–compression test. The behavior of materials under cyclic loading conditions is an issue which might be considered both in micro- and macroscale.

Most of metals, which undergo the plastic deformation, have ability to strain hardening. The hardening is defined as an increase in the yield stress caused by an increase in the plastic deformation [5]. The hardening effect results from the dependence between the yield stress and the history of the plastic straining [6]. In numerical simulations of the material hardening behavior, different hardening models might be used. The isotropic models assume a uniform expansion of the yield surface without its translation in the stress space. The translation of the yield surface without the change in its size is described by kinematic hardening models. The Frederick–Armstrong (F-A) model [7] is a commonly used nonlinear kinematic hardening model which includes the Bauschinger effect. In the Chaboche kinematic hardening model, the back stress is the superposition of several kinematic laws [8]. The Chaboche model enables a better description of the material response to the external load than F-A model and is preferred in simulations of cyclic load problems, e.g., shakedown or ratcheting [6]. Other popular kinematic hardening models are Ohno–Wang, Teodosiu–H and Yoshida–Uemori [9]. The combination of nonlinear kinematic and isotropic hardening laws was included in Chaboche and Lemaitre model [10].

The determination of hardening parameters and modeling of hardening curves are very important problems under scientific consideration. In [11], the usefulness of the finite element model updated procedure (FEMU) for the identification of hardening parameters of the 2024-T3 aluminum alloy was examined. The determination of hardening parameters by means of the virtual field’s method (VFM) was tested in [12]. The influence of the number of tension–compression cycles on the values of hardening parameters with the use of Chaboche and Lemaitre mixed hardening model was included in [6].

The determination of material parameters can be also based on the micro–macro identifications. Placidi et al. [13] have identified the set of parameters which characterized the model of homogeneous linear second gradient D4 orthotropic continuum. In this research, the set of constitutive parameters includes: the material of fibers, fibers cross sections and the internal rotational springs positioned at each intersection between the two groups of fibers. Misra and Poorsolhjouy [14] have used the continuum model of the granular media enriched by non-classical terms in order to describe the material behavior more precisely. The higher-order elastic constants of the enhanced continuum model were determined from discrete simulations subjected to relevant boundary conditions. After that, the elastic constants were used for the determination of the microscale stiffness coefficient of the continuum model. In [15], the Dahl and LuGre considered models which describe the friction in the concrete specimens subjected to uniaxial compressive loading. De Angelo et al. [16] showed the usefulness of homogenized models for the description of the pantographic microstructure of metallic specimens. Detailed information concerning the aforementioned research is included in [13,14,15,16].

The proper selection of hardening parameters plays an important role both in experimental and numerical calculations. However, material hardening parameters obtained from experimental investigations and/or from numerical simulations can be found with a certain precision. Several samples made from the same material can provide little different hardening parameters. Nowadays, soft-computing methods, e.g., neural networks and genetic algorithms, become popular approaches which allow improving the magnitudes of material hardening data. In this paper, a general approach to enhance material parameters by means of the fuzzy set theory is proposed.

The advantage of the fuzzy set theory is that it can be applied in the determination of hardening parameters assuming their scattering. The fuzzy logic is a very efficient non-deterministic computational method for solving engineering problems for which the knowledge about considered model is incomplete or even unknown. The fuzzy theory and fuzzy sets are described in [17]. In [18], the application of the fuzzy logic for solving simple mechanical problems was presented. The response of the Bodner–Partom material to external loads with the use of fuzzy logic was investigated in [19]. The application of the fuzzy set theory for the determination of material strain hardening parameters in Hollomon equation based on the Heyer’s test was described in [20]. In [21], the fuzzy logic model was also used for the prediction of failure of composite plates under different temperatures. Akkurt et al. [22] applied the fuzzy logic model for the prediction of the compressive strength of cement mortar. The fuzzy set theory has been also applied in commercial applications. Vygnanov and Fironov [23] proposed the fuzzy logic system estimating the risks of the high-speed railway. Guclu and Metin [24] tested the fuzzy logic system to control the light rail transport vehicles vibrations. The application of the fuzzy logic for the durability assessment of bridge structures is presented in [25].

In this paper, the application of the fuzzy set theory for adjusting material hardening parameters found from experimental stress–strain curves obtained in the cyclic tension–compression tests is presented. The \(\upalpha \)-level optimization method was used as the main numerical tool. Material hardening parameters are assumed to be fuzzy variables. Numerical hysteresis loops are compared to experimental ones, and the fuzzed error of approximation was found. In the defuzzification process, the most reliable hardening parameters are found. The advantage of using the fuzzy set theory over typical statistical methods is that the fuzzy logic considers not only the scattering of input and output parameters but also the usually nonlinear properties of the mapping model (the plastic flow rule). The future goal of this research is the analysis of sophisticated elastic–plastic material models like generalized elasticity and plasticity models which includes more nonlinear effects, e.g., creep and relaxation, as well as, visco-plastic models, etc. In the further research concerning the modeling of the hardening curves, the Bodner–Partom model (generalized plasticity) as well as the Chaboche model will be tested.

Before such a study, the classical elastic–plastic nonlinear models are considered, namely Voice isotropic and Frederick–Armstrong kinematic models. It might be treated as a first stage of a research concerning the usefulness of fuzzy logic for the determination of hardening parameters.

2 Combined hardening: problem formulation

In the numerical simulation of elastic–plastic stress–strain curves, the Hooke’s law and the plastic flow rule associated with von Mises yield criterion were used. For simplicity, the explicit integration was used. The basic equations for combined isotropic–kinematic hardening model are presented below [26].

The yield criterion is described by the following Eq. (1) [27]:

$$\begin{aligned} \hbox {f}(\varvec{\sigma }, \mathbf{x}, p) = \sqrt{\frac{3}{2}(\varvec{\sigma }'-\mathbf{x}) : (\varvec{\sigma }' - \mathbf{x})} -\sigma _y -r (p) \end{aligned}$$
(1)

where \(\varvec{\sigma }\) and \(\varvec{\sigma }'\) are the stress and deviatoric stress, respectively, \(\mathbf{x}\) is the back stress, \(\sigma _y \) is the yield stress, r(p) is the isotropic hardening function, and p is the effective plastic strain. Elastic deformation occurs if the yield function \(\hbox {f}(\varvec{\sigma }, \mathbf{x}, p)< 0\), while the plastic flow takes place for \(\hbox {f}=0\) [28, 29].

The effective plastic strain is described by the following formula (2) [30]:

$$\begin{aligned} \mathrm{d}p= \sqrt{\frac{2}{3}\mathrm{d}\varvec{\varepsilon }^{\mathrm{p}} : \mathrm{d}\varvec{\varepsilon }^{\mathrm{p}}} \end{aligned}$$
(2)

where \(\mathrm{d}\varvec{\varepsilon }^{\mathrm{p}}\) is the increase in the plastic strain.

The consistency condition is described by Eq. 3:

$$\begin{aligned} \mathrm{d}f= \frac{\partial f}{\partial \varvec{\sigma }}: \mathrm{d}\varvec{\sigma }+\frac{\partial f}{\partial \mathbf{x}} :d\mathbf{x}+\frac{\partial f}{\partial p}\mathrm{d}p=0 \end{aligned}$$
(3)

The normality condition is (4):

$$\begin{aligned} \mathrm{d}\varvec{\varepsilon }^{\mathrm{p}} = \mathrm{d}\lambda \frac{\partial f}{\partial \varvec{\sigma }} \end{aligned}$$
(4)

where \(\mathrm{d}\lambda \) is the plastic multiplier which can be found from the consistency condition (3). The plastic flow associated with the von Mises yield criterion is described by Eq. 5 [27]:

$$\begin{aligned} \mathrm{d}\varvec{\varepsilon }^{\mathrm{p}} = \mathrm{d}\lambda \frac{3\varvec{\sigma }'}{2\sigma _e} \end{aligned}$$
(5)

where \(\sigma _e \) is the effective stress. Otherwise, \(\mathrm{d}\varvec{\varepsilon }^{\mathrm{p}}\) equals 0 [28].

The back stress increment \(d\mathbf{x}\) in Frederick–Armstrong nonlinear kinematic hardening model [31] is:

$$\begin{aligned} \mathrm{d}{} \mathbf{x}=\frac{2}{3}c\mathrm{d}\varvec{\varepsilon }^{\mathbf{p}}-\gamma \mathbf{x}\mathrm{d}p \end{aligned}$$
(6)

where c is the material constant, \(\gamma \) defines the rate of the kinematic hardening saturation, and \(\mathrm{d}p\) is the equivalent plastic strain increment.

The increment in isotropic hardening dr(p) is determined by Voce isotropic law [32] as:

$$\begin{aligned} \mathrm{d}r(p)=b(Q-r)\mathrm{d}p \end{aligned}$$
(7)

where Q is the saturated value of r and b determines the rate at which the saturation is achieved.

Assuming the superposition of elastic and plastic strain increments \(d \varvec{\varepsilon } =d \varvec{\varepsilon }^{\mathrm{e}} +d \varvec{\varepsilon }^{\mathrm{p}}\), the Hooke’s law is described by the following Eq. (8):

$$\begin{aligned} \mathrm{d} \varvec{\sigma } =\mathbf{C} (d \varvec{\varepsilon } -\mathrm{d}\varvec{\varepsilon }^{\mathrm{p}}) \end{aligned}$$
(8)

In (8), the Voigt notation is used for stress and strain tensors and \(\mathbf{C}\) is the constitutive matrix. Combining Eqs. (1-8), one can express the plastic multiplier in terms of the strain increment:

$$\begin{aligned} \mathrm{d}\lambda =\frac{\mathbf{n}\cdot \mathbf{C}\mathrm{d}\varvec{\varepsilon }}{\mathbf{n} \cdot \mathbf{Cn}-\gamma \mathbf{n}\cdot \mathbf{x} +\frac{2}{3}\hbox {c}{} \mathbf{n}\cdot \mathbf{n}+b(Q-r)} \end{aligned}$$
(9)

where \(\mathbf{n}=\frac{\partial f}{\partial \varvec{\sigma }}\). It should be mentioned that for von Mises yield criterion, \(\mathrm{d}\lambda =\mathrm{d}p\).

The plastic multiplier can be also expressed in terms of the stress increment. After substituting (8) into (9) and some transformations, the plastic multiplier is described by Eq. 10.

$$\begin{aligned} \mathrm{d}\lambda =\frac{\mathbf{n}\cdot \mathrm{d}\varvec{\sigma }}{\frac{2}{3}c\mathbf{n}\cdot \mathbf{n} -\gamma n\cdot \mathbf{x} +b(Q-r)} \end{aligned}$$
(10)

Equations (1) to (8) were used in developing of numerical algorithm for an elastic–plastic material response to the kinematic load. For assumed explicit integration, the algorithm executes as follows:

  1. 1.

    For a given strain increment, check the yield condition (1)

  2. 2.

    If \(f < 0\), find the stress increment from (8)

  3. 3.

    If \(f=0\),

    • find the plastic multiplier from (9)

    • find the plastic strain increment from (5), the back stress increment from (6) and the isotropic hardening increment from (7)

    • find the stress increment from (8)

    • upgrade strain, stress, back stress and isotropic hardening at the end of this step.

For a simple uniaxial loading, the explicit form of the stress increment vs. the strain increment for \(f =0\) can be derived (11):

$$\begin{aligned} \mathrm{d}\sigma =E\left( 1-\frac{E}{E+c-\gamma \cdot x+b(Q-r)}\right) \mathrm{d}\varepsilon \end{aligned}$$
(11)

where E is the Young’s modulus.

3 General description of the fuzzy logic analysis

Fuzzy logic allows including various uncertainty of the input and output data [33] in an analysis. In general, it represents a mapping of fuzzy input variables into the result space. The determination of the fuzzy output variable includes some steps among which the most important are the fuzzification and defuzzification (Fig. 1). In the fuzzification step, the membership functions associated with input variables should be defined. After that, the fuzzy output variables are determined with the application of the mapping model [34]. The mapping model might be determined explicitly by the function, can be described by means of linguistic variables or can be a sophisticated algorithm, e.g., the system of differential equations.

Fig. 1
figure 1

Steps of the fuzzy analysis

Full size image

A number of fuzzy logic methods can be applied in order to obtain the membership functions for fuzzy output variables. In this research, the \(\upalpha \)-level optimization as well as the extension principle method was tested. In the extension principle, the membership function for z output variable is calculated with the use of minimum and supreme operators according to the following formula (12):

$$\begin{aligned} \mu (z)=\sup \min [\mu (x), \mu (y)] \end{aligned}$$
(12)

where \(\mu (x)\) and \(\mu (y)\) are membership functions of x and y input variables, z is an output variable, and \(\mu (z)\) is the membership function of z output variable.

For randomly selected x and y fuzzy input variables, their membership functions \(\mu (x)\) and \(\mu (x)\) are computed. After that, the minimum values of \(\mu (x)\) and \(\mu (x)\) are determined. In this way, z variable can be associated with several \(\mu (z)\) values. Thus, the supremum operator is used. For floating-point numbers, the maximum \(\mu ({z})\) value should be searched for in the \(\hbox {z}\pm \Delta \hbox {z}\) interval.

The main disadvantage of the extension principle method is the favoring of minimum values of \(\mu (z)\). (The minimum operator is applied at first.) Relatively smooth envelope (an application of supremum operator) can be obtained only for a very dense searching for x and y fuzzy input variables. This is time-consuming process. The results obtained by means of the extension principle usually receive a form of a jagged solution. It is worth noting that the final result also depends on the assumed precision \(\Delta z\).

If the mapping function requires a lot of calculations, the extension principle is not effective. Much better results might be obtained by the application of \(\upalpha \)-level optimization method [17].

Fig. 2
figure 2

The \(\upalpha \)-level optimization method

Full size image

In this approach [18], the membership functions of x and y fuzzy input variables are separated into sufficient number of \(\upalpha \)-levels. For the specific \(\alpha _k \) level \(\alpha _k\) [0; 1] the boundary values \(x_{\alpha kl}\), \(x_{\alpha kr} \), \(y_{\alpha kl} \) and \(y_{\alpha kr}\) are computed (Fig. 2). For randomly selected \(x\in [x_{\propto kl}, x_{\propto kr}]\) and \(y\in [y_{\propto kl} ,y_{\propto kr}]\), the minimum \(z_{\alpha kl}\) and maximum \(z_{\alpha kr}\) are found. The above procedure is repeated for all \(\alpha \)-levels. Extreme z values for all \(\alpha \)-levels determine the shape of the fuzzy output membership function \(\mu (z)\).

The \(\upalpha \)-optimization method can be used only if the mapping operator is continuous and the fuzzy result space is convex. In comparison with the extension principle, the \(\upalpha \)-level optimization method enables to obtain a smooth shape of the output membership function with the same number of calculations. For this reason, the \(\upalpha \)-optimization method was used as the main numerical tool in this research. The extension principle was applied only in order to check the consistency of the obtained results.

In the defuzzification procedure, the most valuable result is selected from the fuzzy result space. There are several defuzzification methods, e.g., the height method, the mass center method, level rank method or Jain and Chen method. In this research, the mass center method was applied. This concept is based on the searching for the center of the space below the fuzzy output membership function. Mass center is defined by the following formula [17]:

$$\begin{aligned} z_0 =\frac{\int z\cdot \mu (z)\mathrm{d}z}{\int \mu (z)\mathrm{d}z} \end{aligned}$$
(13)

where \(z_0\) is the most reliable value of the z variable and \(\mu (z)\) is the membership function of z variable.

4 Experimental investigations

In the experimental investigations, the hardening curves for samples subjected to cyclic tension–compression loads are recorded. The cyclic tests with certain strain amplitudes enable to identify material hardening parameters [35].

The strain-controlled cyclic tension–compression test was carried out at room temperature on ZWICK/ROELL Z030 testing machine. The deformation of \(\pm 5\%\) of measuring base was applied. The sample was made of PA6 aluminum alloy. The specimen geometry is shown in Fig. 3.

Fig. 3
figure 3

Dimensions (mm) of the specimen used in this experimental research

Full size image

The experimental curve for PA6 alluminium alloy subjected to the tension-compression test is presented in Fig. 4. Engineering quantities \(\sigma _\mathrm{nom}\) and \(\varepsilon \) registered in the test are transferred into true ones by formulas (14) and (15) assuming material constant volume for the plastic deformation.

$$\begin{aligned} \sigma= & {} \sigma _\mathrm{nom} (1+\varepsilon ) \end{aligned}$$
(14)
$$\begin{aligned} \phi= & {} \hbox {ln}(1+\varepsilon ) \end{aligned}$$
(15)

where \(\sigma _\mathrm{nom} \) is the engineering (nominal) stress, \(\varepsilon \) is the engineering strain, \(\sigma \) is the true stress, and \(\phi \) is the logarithmic strain.

The true stress in (14) is defined as the force related to the current cross section. The logarithmic strain (15) enables the superposition of strains in the large deformation analysis. A material subjected to finite deformations will have the same stress–strain behavior in tension and compression if the true stress and the logarithmic strain are applied.

Fig. 4
figure 4

The experimental curve for PA6 aluminum alloy

Full size image

5 Numerical simulations of the hysteresis loops

The equations presented in Sect. 2 are used in order to develop the numerical program which computes the elastic–plastic material response to the assumed strain increment. This program implemented in MATLAB simulates the cyclic tension–compression test. For simplicity, the explicit integration procedure is applied with about four thousand integration steps for single hysteresis loop. (Explicit procedure is conditionally stable.) For initially selected \(b, Q, c, \gamma \) hardening parameters, the stress–strain curve is found as shown in Fig. 5.

Fig. 5
figure 5

Comparison of experimental (a) and numerical (b) curves for a cyclic tension–compression test

Full size image

One can see from Fig. 5 that quite a good agreement between experiment and numerical simulation was achieved. In the initial selection of hardening parameters, the following indicators are used:

  • the size of the initial yield surface \(\pm \sigma _{Y0}\) and the size of the yield surface for the final loop—determination of Q parameter

  • a number of loops up to hysteresis stabilization—parameter b

  • extrapolated shift of the yield stress for the final stable loop—the ratio \(\frac{c}{\gamma }\)

  • curvature between elastic and plastic response—parameter \(\gamma \)

As the initial guess, the following hardening parameters are selected: \(Q = 154\ \hbox {MPa}\), \(b =7\), \(c = 12{,}000\ \hbox {MPa}\), \(\gamma = 100\) (\(\frac{c}{\gamma } = 120\ \hbox {MPa}\)).

Better fitting of numerically simulated curve to the experimental one can be achieved by means of the optimization approach [36]. Here, the following error norm (B) was introduced:

$$\begin{aligned} \Vert B\Vert = \sqrt{\int _{\varepsilon } (\sigma _\mathrm{exp} - \sigma _\mathrm{app})^{2}\mathrm{d}\varepsilon } \end{aligned}$$
(16)

where \(\sigma _\mathrm{exp}\) and \(\sigma _\mathrm{app}\) are experimental and approximation stress values, respectively.

In the optimization process, all hardening parameters are randomly distorted up to 20% and for each set of \(b, Q, c, \gamma \) parameters, the error norm (15) is computed. The minimum error (\(B=\) 12.65%) is reached for the following data: \(Q =138.42\ \hbox {MPa}\), \(b = 7.82\), \(c = 9617.90\ \hbox {MPa}\), \(\gamma = 118.91\). The comparison of numerical and experimental hysteresis loops is shown in Fig. 6.

Fig. 6
figure 6

Comparison of experimental and numerical curves after the optimization of hardening parameters

Full size image

It can be clearly seen in Fig. 6 that very good fitting of both curves is reached, even for only four hardening parameters taken into consideration. Better fitting requires using more complex hardening models, such as Chaboche and Lemaitre model.

6 The application of fuzzy logic in enhancement of hardening parameters

The optimization procedure mentioned in the previous section provides a very good agreement between cyclic tension–compression experiment and its numerical simulation. However, similar (small) error norms (16) can be obtained for various sets of hardening parameters. The best fit can be applied for the first loop and for the last loop, or a mild fit can be applied to all loops simultaneously without the preference of the particular one. In all these cases, the integral will be similar. Therefore, the question arises, which set provides the best simulation of the material response to the external load.

Fuzzy logic analysis is the tool which can help to select the best set of material hardening parameters assuming some uncertainty of the hardening data. In this research, the hardening parameters for combined hardening model (Q, b, c, \(\gamma \)) are assumed to be fuzzy input variables. These parameters are randomly distorted up to \(20\%\). The triangular membership functions for each set of the input parameters are built. After that, the fuzzy analysis is performed. The fuzzy output variable \(z\equiv B\) and its membership function \(\mu (z)\) are calculated. The most reliable \(z_0\) value (crisp value) is found in the defuzzification process. It is assumed that the best set of material hardening parameters provides the minimum of discrete \(z_0\) value.

This procedure is time-consuming and may take several hours of computation depending on the number of the data randomizations, number of \(\upalpha \)-levels and number of integration steps for a single hysteresis loop.

As the result of the fuzzy logic analysis, the following hardening parameters are computed: \(Q =134.31\ \hbox {MPa}\), \(b =5.98\), \(c =9925.36\ \hbox {MPa}\), \(\gamma =101.22\). The approximation error B calculated by Eq. 16 is 27.96%. The comparison of numerical and experimental hysteresis loops is shown in Fig. 7.

The advantage of using fuzzy logic over classical methods of statistics is that the fuzzy set theory considers the randomness of both input and output data. The nonlinear behavior of the mapping procedure is also included in the analysis.

Fig. 7
figure 7

Comparison of experimental and numerical curves after the fuzzy logic analysis

Full size image

7 Summary and conclusions

The proper selection of hardening parameters is necessary both for experimental and numerical calculations. In this paper, the procedure for the selection of hardening parameters of mixed nonlinear kinematic and isotropic hardening model with the use of authorial algorithm was implemented. The experimental as well as numerical hysteresis loops, for which hardening parameters were selected, were compared with the use of least-squares method. The analysis of results showed a relatively good fitting of both curves.

In order to improve the results of numerical simulation, the fuzzy set theory was applied. Fuzzy logic is a non-deterministic method for the analysis of data which is imprecisely defined. The application of fuzzy logic in plasticity theory is associated with the fact that many parameters are usually imprecise or uncertain. An example of this would be load magnitudes, material data and the friction coefficient which might change from one product to another. Therefore, the application of fuzzy logic for solving engineering problems for which input parameters are imprecisely defined is deliberate.

Nonlinear hardening models are good examples of fuzziness. The hardening parameters might be determined in an experimental test and optimized with the use of a numerical procedure. Better results might be achieved by the execution of the fuzzy analysis. In this research, the hardening parameters of mixed nonlinear kinematic and isotropic models were selected with the use of fuzzy logic. The variation of hardening parameters was assumed, and the fuzzification was done. As a mapping model, a system of differential equation associated with plasticity theory was used. Many numerical simulations were carried out, and the lowest discrete value obtained in a defuzzification step was selected as a final result.

The simulation shows the potential of application of fuzzy logic analysis for the determination of hardening parameters. The fuzzy set theory allows finding the variation of elastic–plastic material response when uncertainness of input parameters is taken into consideration. Therefore, the fuzzy analysis might be a useful tool in the analysis of metal forming processes.