FEM Simulation of Non-proportional Multiaxial Fatigue Damage

. The paper deals with implementation of the optimizing process into multi-axial rainflow analysis and cumulative damage calcula tion. It’s presented computational program FEA_FAT compiled in MATLAB. Stress analysis is realized by finite element procedure and the cumulative damage can be calculated by using two fundamental ways – critical plane approach and so called integral approach. Testing example presents random stress analysis and damage prediction of a simple FE model with non-proportional loading.


Introduction
The degree of stress or damage of structural elements is most often determined by evaluating the stress or strain response. These are regularly described by time-dependent stress tensor (relative strain), i.e. not by one but by several components. This leads to the question -how to determine the accumulation of damage in this case. If the placing, geometry and load of the object show some degree of proportionality (Fig. 1), it is possible to accept modified uniaxial "quasi-static" strength hypotheses (Tresca, von Mises) to calculate the equivalent stress (relevant for the damage accumulation calculation). This assumption is not particularly satisfactory in cases of random and often resulting disproportionate (Fig. 2) and phase shifted loading. Then uniaxial hypotheses may not suit reality. Especially in the case of finite element analyzes we get a stress tensor, which in most tasks contains more than one significant component. Then there is a problem how to realize rainflow decomposition for the whole stress tensor. The application of a classical equivalent stress is not directly possible, since we only get positive values and counting closed cycles would not be justified [7,9,11,12]. Other procedures must be designed to meet the real state of fatigue damage estimation of the investigated place. Several procedures were proposed for multiaxial rainflow decomposition and subsequent estimation of damage rate in the past. We have chosen and analyzed two basic approaches in our study -the socalled. Integral Approach (IA) and Critical Plane Approach (CPA). There are many , different variants and attempts in the literature to define criteria for determining the degree of fatigue damage. The best-known criteria are the Findley and Dang Van criterion of the critical load plane, the Sines and Crossland criterion based on stress invariants, the Papadopoulos criterion of stress average in elementary volume, the Garund criterion based on energy concepts, the integral method of linear combination of stress components and others [1][2][3]5,6].

Design and Analysis of Methodologies for Fatigue Damage Determination by FEA Application
Based on the analysis of the excitation of the investigated object, we can build a static or dynamic computational finite-element model (we consider linear now). In the case of analysis of static or quasi-static time-varying loading, we can solve the problem directly by solving equilibrium equations based on FEM principles and by calculating the necessary stress response. If only the amplitude of the load changes, then it is possible to use the analysis for the so-called unit loading with subsequent calculation of unit stress. Its real values are obtained by superposition from contributions of individual loading amplitudes. We apply the known equation of motion if the inertia and damping effects play an important role. By its solution we get the displacement field and consequently the stress field. This process can again be accomplished by directly integrating the equations of motion using any of the known numerical methods. It is a universal method that does not depend on the degree of non-linearity of state equations. Another approach is based on the principle of the above-mentioned unit load superposition. The advantage is that if the load amplitude is changed, we get the necessary outputs quickly. The disadvantage is the applicability only for linear or some linearized tasks. When applying FEM, an approach based on the known modal decomposition of the state equation and subsequent calculation of the so-called modal stresses is often used. From these, real values can be obtained relatively quickly. There are several options and it is up to the analyst's experience to choose from them. As already mentioned, the biggest problem in the case of several significant components of the stress or strain tensor remains the need for further processing of the equivalent stress value (or deformation) -most often by the rain flow method. There are many approaches based on the proportionality of the load process. They give acceptable results to the extent possible [1,10,14]. However, this is worse in the case of significant disproportionality in the load, where their applicability is limited or only indicative. The paper deals with selected methods of critical plane approach (CPA) as well as with two modifications of the integral approach, which were implemented in FE analyzes [13,15,19,21].

Findley's Method of Critical Plane
Findley assumed that the critical plane in terms of fatigue crack initiation would be the plane with the greatest shear stress [7,12]. He defined the equivalent shear stress in the form where k is the coefficient that the author has determined for the tough material to be around 0.3. If we consider that we have only a 1-axis alternating tensile-pressure test, we will adjust the previous relationship with a certain amount of von Mises simplification as follows ( The equivalent stress prepared in this way can be distributed by rainflow decomposition. The approach is particularly suitable for proportional loads, but it is also applicable for nonproportional loads. A certain problem is the appropriate choice of the coefficient k. For example, in [15] it is stated that Findley's suggested value k = 0.3 is rather exaggerated, which is also documented by our numerical tests.

Dang Van's Critical Plane Method
Dang Van also assumed that the critical plane for fatigue crack initiation would be the plane with the greatest shear stress. He defined the equivalent shear stress as where C is the fatigue limit for torsion, C is the fatigue limit for axial tension-pressure and 1, 2, 3 are the principal stresses. Consider again that we only have a 1-axis alternate tension-pressure test, and then we can use von Mises to adjust the previous relationship to In relation (4) we can still make a simplification adjustment: , Prepared equivalent stress (4) can already be decomposed into closed cycles by rainflow decomposition. This approach is more appropriate for the proportional load [17,20]. However, its use is acceptable at a lower level of disproportionality based on the experience of the authors or gives the possibility to orientate in the values of damage accumulation.

Numerical Method of Critical Plane Search for Modified von Mises Stress
Another method is searching the critical plane in a numerical way, i.e. finding an nCPA normal vector of the plane that is critical to fatigue crack initiation. The equivalent stress is defined as follows where  (nCPA) is normal and  (nCPA) is the shear stress in the critical plane determined by the nCPA vector. To find nCPA, one of the well-known MATLAB search procedures known as optimization methods can be used. This means that we can use procedures such as fmins.m, fminsearch.m, patternsearch.m and more.

Integral Method of Linear Combination of All Stress Components
We consider the n-dimensional vector of stresses  (t) = [x, y, ...] T with a random character as a local response to the load by common and random external forces. The basic idea of the method is to define the equivalent stress as a linear function of the individual components of the tensor, i.e. where . This defines a spherical surface in the coordinate system c1, c2, c3. The numerical task is to find on this surface a combination of the components of c vector for which the accumulation of damage will be greatest. In the case of spatial stress, c vector contains six components. This can complicate numerical analysis over time. The search for a suitable combination of c vector components means solving an extremelization task linked to a six dimensional hypersurface. Therefore, other techniques have been proposed to minimize the dimension of the c vector. It is necessary to recall the appropriateness of using the presented approach also in cases of non-proportional loading.

Integral Method of Linear Combination of Principal Stresses
As mentioned above, the search process efficiency of the c vector components is complicated by their number. One solution for reducing the number of components ci is to perform all calculations for the principal stresses. This also leads to the problem of finding three unknown coefficients of ci in the case of spatial stress. The equivalent stress in terms of multiaxial rainflow decomposition will be calculated as follows where 1  2  3 are principal stresses. This procedure is well applied to a nonproportional loading too. Based on the authors' experience, it should be noted that IA with principal stresses (3-dimensional task) underestimates the value of damage accumulation over the previous 6-dimensional approach. When used for proportional load cases, the calculated fatigue damage rate is less than for CPA methods.

Integral Approach Implementation to the Fatigue Analysis by FEM Application
Similarly to the CPA method of finding the critical plane for modified von Mises stress, the IA methods are used in the FE analysis of virtual models too. Authors therefore consider the random excitation of the investigated structure, the geometrically and materially linear computational model and the Miner's hypothesis of the linear accumulation of fatigue damage. Under these assumptions, we will present an effective numerical procedure for quickly finding of c vector. The basic idea of numerical search for the "critical" values of c vector is to define an optimization task to find maximum damage D. According to the authors, this is the fastest way to find such a combination of ci for which D  max. The goal will be to find optimal c'vector values for the target function in the form where t is the implementation time interval, D is the degree of damage accumulation, Ni is the number of cycles to damage at the i th stress level determined from the S/N curve, ni is the number of load cycles at the i th stress level. Using rainflow decomposition of equivalent stress MRF (t, c') we get ni and we can also apply one of the hypotheses of damage accumulation calculation. Based on the theoretical considerations, the FEA_FAT calculation program was created in the MATLAB programming language [16,18]. The program combined stress FEM analysis and node damage accumulation calculations based on maximizing D damage. Methodical design of multiaxial calculation of fatigue damage (or fatigue life) is presented in Fig. 3. Brief characteristics of individual programs and calculation procedures of the developed software package:  FEA_FATmain control program -it reads the set of input data, processes the excitation vector for the whole time interval, after processing the stress behavior in the elements it performs damage calculation and predicts fatigue life for each element by means of the UNAVA_NODOPT program,  STATICsolver of static or quasi-static tasks, calculates displacements and stresses in elements or nodes of the model,  DYNAsolver of dynamic tasks solved in time domain by application of Crank-Nicolson numerical method, calculates time dependence of displacements and stresses in elements or nodes of the model,  STOCHASTsolver of stochastic oscillation in the frequency domain, it calculates the power spectral density (PSD) as well as the covariance matrix of displacements and stresses in the elements or nodes of the model,  UNAVA_NODOPT, MCHRF, DAMAGE, FMINSEARCHsubroutines that processed both 1-axis and multi-axis fatigue damage problems using the proposed optimization approach for multiaxial rainflow decomposition and calculate the damage rate and estimated time to damage,  DAMAGE_Rsubroutine that calculates the estimated time to damage from the PSD stress in each element by Rajcher's methodology. The proposed methodology includes an analytical-computational view of the process of loading, modeling and method selection in the task of predicting the accumulation of fatigue damage -especially in cases of multiaxial stress [4,8,11,13]. The procedures and obtained information correspond to the results of the work of other authors. Finding a universal or best hypothesis for fatigue is a problem that is likely to be answered by future generations of experts.  (Fig. 4). The chosen results are presented in Table 1 and 2.

Conclusions
The paper presents an original computational approach of fatigue damage accumulation prediction based on simulation analysis using special search algorithms built in MATLAB and on principles of FEM. Selected known and less known methodologies were presented and compared. Three were based on the principle of finding a critical plane from the perspective of fatigue crack initiation and the fourth (IA) was based on a linear combination of stress tensor components. Although the physical nature of IA is somewhat vague and based rather on a purely mathematical description of the problem, IA is gaining more and more supporters in the professional community. The main aim of the authors was to design and present the original methodology and algorithm (implemented in MATLAB) for the application of optimization algorithms in the computationally demanding process of fatigue damage prediction. All this was done by means of a finite element analysis assuming the random and non-proportional character of the stress response. The obtained results confirm the authors' reflections on a suitable alternative to more popular methods based on the search for the critical plane.