EFFECTIVE ALGORITHM FOR STRUCTURAL OPTIMIZATION SUBJECTED TO FATIGUE DAMAGE AND RANDOM EXCITATION Summary

The aim of the paper is to present the computationally timeconsuming task of simulating the process of randomly oscillating thin-shell structures to realize an optimal design with limitations in terms of fatigue damage. The most important aim here is to design an effective optimization algorithm and choose an appropriate approach for the evaluation of multiaxial fatigue damage in the context of the random and non-proportional character of stress in the structure when considering the standard finite element model. The authors propose their own optimization algorithm, which is described in the present study and, on the basis of tests, has proven to be suitable for the aforementioned technical problems. The proposition of algorithms for calculating the accumulation of fatigue damage in non-proportional and multiaxial stresses (especially in terms of the application of rainflow analysis by decomposition of the equivalent stress, as determined by the appropriate “fatigue” criterion) is very important in such


INTRODUCTION
It is almost impossible to pick up a journal or conference programme focused on computational mechanics that does not contain some reference to structural optimization.Although designing machine parameters by experience is achievable, it is much better and more effective to predict the base properties of the new designed structure by using an optimization procedure, which is generally predicated on a series of controlled computing analyses [5].
Today, we expect designed objects to be optimally balanced over the entire life cycle, i.e., projecting, manufacturing, running, maintaining and liquidating.The mentioned processes mainly relate to the economic aspects of each stage.Achieving this goal is very difficult, because designers are usually confronted by contradictory demands related to the individual stages of the aforementioned life cycle of designed objects [3,8,11].
The first formulations of optimization problems in the context of mathematical programming appeared from around 1960.One of the pioneers who significantly influenced the development of optimal construction designs for machines and their components was undoubtedly Schmit.He linked optimization methods with the new and progressive computational methods of the time, which included the finite element method.During this period, the weight of the monitored object or a certain strength condition was the objective function.Optimization processes were gradually improved by adding other limiting conditions.In the second half of the last century, other works of a similar nature, which extended options in the field of optimally designing the parameters of machines and their components towards automated approaches, were published.We cannot omit the works of Kirch, Morrow or Gallagher here.Numerous effective approaches were proposed, based not only on a purely mathematical comprehension of the optimization problem, but also on slightly non-traditional or more precisely unaccustomed approaches, which play an important role in solving various technical problems [9].These approaches use some of the basic principles of mechanics.
The optimization of mechanical systems combines numerical mathematics and engineering mechanics.It is used in applications in civil engineering, mechanical engineering, the automotive and shipbuilding industry, etc.It has made the biggest progress in the last 30 years thanks to the utilization of very fast numerical computers and computer graphics.When choosing cost, structural weight or maximum power at a limited cost as the design criterion, the importance of optimization is evident [9].

STRESS ANALYSIS USING THE FINITE ELEMENT MODEL
Generally, the stress calculation for the selected element or node of the finite element model can be realized as follows [2,7]: 151.
-A) In the case of the quasi-static time-dependent load, we must solve the following equation: (1) while, for the j-th element (node), it is possible to calculate the stress response using the relationship: where Dj is the standard material matrix, Bj is the transformation matrix between the strain function and the j-th element node displacements, Tj is a classic transformation matrix between the local and global coordinate systems, and Aj is a Boolean matrix, i.e., the localization matrix determining the element displacements' position in the global displacement vector.
-B) In the case of a linear oscillating system described by the equation below: ) it is possible to find a solution by applying a well-known modal transformation, (4) of the original Eq. ( 3) to a significantly smaller set of differential equation in the form: (5) where y(t) is the vector of the so-called modal coordinates, V is the modal matrix of the preferred modal shapes (eigenvalues vectors), λ is the corresponding spectral matrix, I is identity matrix, and B ~ is the damping matrix for the applied modal shapes [2,7].Solving Eq. ( 5) provides the time response y (t), which enters the resulting relationship for the stress response in the j-th element as follows [2,7]: ) Previous theoretical considerations were of a general nature.Let us now focus our attention on the well-known shell finite elements (Kirchhoff's or Mindlin's formulation) [2].The stiffness parameters depend on material constants and element geometry (mainly its thickness).At first, we have to prepare the stress calculation process.This process is based on the expression of the j-th element membrane forces and bending moments (without shear forces) [2], i.e.: The integration matrices Im and Ib are: and can only be calculated by using the numerical approach.Vector j L u is the displacement vector of the j-th element in the local coordinate system and t is the element thickness.Further details about the material and element matrices Em, Eb, D, Bm, Bb are presented in the relevant literature [2].The extreme stress values can be expected at the top or on the bottom surface.Generally, this means:  (10) or: Let us build new material and auxiliary matrices as follows: where the matrix I3 is the classical unit matrix.Then (11) can be rewritten as follows: By assuming a relation between the local element displacements and the global displacement vector u(t): (15) (13) and ( 14) may be rewritten as: where TLG is a classic transformation matrix between the local and global coordinate systems, and T01 is again a Boolean matrix, i.e., the localization matrix determining the element position in the global displacement vector u (t).The response vector u (t) can be obtained from the solution of Eq. ( 1) or (3).

Multilevel S-N curve
It is well-known that the Wöhler curve (Fig. 1) or S-2N curve is the basic source of information on material fatigue life.Generally, the S-2N curve is statistically evaluated by an experimental fatigue curve [1,11], which reflects the behaviour of the magnitude of a cyclical nominal stress Sa (or A in subsequent analysis) versus the logarithmic scale of cycles to failure 2Nf.It is advantageous to show it in terms of logarithmical or semi-logarithmical coordinates.The A-2Nf relation can be written as follows: where f is the fatigue stress coefficient, 2Nf is the number of cycles to failure, b is the fatigue strength exponent, and A is the stress amplitude to failure.Some researchers have rewritten the relationship in (18) in the following form [1]: where m=-(1/b) and . Considering the mean stress σm in the modified version of the stress amplitude (using Goodman, Soderberg, Geber), Eq. ( 19) can be rewritten as follows [1]: (20) If k=1 and RF=RE (yield stress), Soderberg's model is used; if k=1 and RF=RM (strength limit), Goodman's model is used; and, if k=2 and RF=RM, Geber's model is used.Using the linear Palmgren-Miner law, we can calculate fatigue damage for stress amplitude Aj as follows [1]: (21) From the experimental measurements, we can construct the so-called multilevel fatigue curve.The mathematical formulation of this multilevel A-2Nf curve for i-th part can be described as follows: where the exponent mi can be written as follows: (23) The calculation of the total damage with respect to the mean stress and multilevel fatigue curve is given by:

Chosen multiaxial damage hypothesis
Let us focus on counting the cumulative damage by using multiaxial rainflow decomposition of the stress response.It should be noted [1,4] that the fatigue damage calculation of the machine parts is generally problematic because the results are considerably changed in the principal stresses [4].Using finite element analysis, we can obtain six components of the stress-time function (multiaxial stress), but it is very difficult to obtain an equivalent uniaxial load spectrum due to comparison with the applied computational fatigue curve.In our case, the rainflow analysis for random stresses, known in its classic uniaxial form as the von Mises or Tresca hypothesis, is impossible.This means that the important goal of this part will be to propose some approaches to estimate the high-cycle fatigue damage for multiaxial stresses caused by random vibrations in the analysed structure [1,4].Generally, the principal approaches for the determination of applicable equivalent stress (or fatigue criterion) in the case of rainflow decomposition are the critical plane approach, integral approach, and the combination of linear stress components; our study focuses mainly on the latter methodology.

Damage calculation using the linear stress components combination approach
The fundamental idea is to count the rainflow cycles for all linear combinations of the stress random vector components [6,10].Practically, if the stress state is biaxial (e.g., thinshell finite element), the stress components can be written in the form of a three-dimensional vector =[x, y, xy] T , such that the equivalent stress will be calculated as follows: . In the case of the shell finite element, we can again obtain following relationships: (26) Hence, the next goal will be to find the extreme value of the estimated damage for vector c and i-th element, i.e.: where Di-max/MRF is the maximum value of the cumulative damage for the i-th element, 2Ni is the number of cycles to failure, and mc is the number of cycles after rainflow decomposition of the stress.Naturally, we have to observe the normality condition for c using the following transformation: The searching process is realized by the FAT_MRFA computational program developed in MATLAB by authors of the article.The program calculates the elements' (or nodes') damage from a stress response using the original optimizing multiaxial rainflow procedure suggested by one of the authors of the article.

Formulation of optimization problem
The optimization problem of the mass minimization, which is subjected to the prescribed fatigue damage or life, is topical.For a structure of multiple elements, the optimization problem of discrete variables [5,8] can be stated mathematically as follows: which subjected to: , k=1, 2, ..., m, (30) where n is the number of elements, m is the number of element groups, Dp is the prescribed cumulative damage, TP is the prescribed fatigue life in hours, D k max is the calculated extreme value of the cumulative damage for k-th element group, and T k min is the calculated extreme fatigue life in hours for the k-th element group.Let us form a new penalized objective function thus: where the penalty function can be:

Optimizing algorithm proposal
The penalized objective function   is solved by new discrete optimizing algorithm (see Fig. 3).The iterative optimizing process is as follows: Fig. 3. Computational scheme of the suggested algorithm Definition: 1. numvar -number of optimization variables 2. Dp -cumulative damage limit 3. V -vector of the discrete optimization variables 4. hran -matrix of the limited values of optimization variables 5. d0 -starting vector (serial number of the vector V values) Finite element analysis and cumulative damage counting for vector V (d0)

Checking the convergence condition END
Prediction of the new optimization variables (serial numbers of the vector V members) using the following iterative relation: Checking the limited values and modification of the vector d (j+1) Finite element analysis and cumulative damage counting for vector V (dj+1) Correction of the vector d (j+1) in the case of an overload of the prescribed limited cumulative damage Dp using the following iterative relation: Checking the limited values and modification of the vector d (j+1) Finite element analysis and cumulative damage counting for corrected vector V (dj+1) where i j x is a serial number of the i-th design variable in j-th iteration step, ceil (Y) is a MATLAB function, which rounds the elements of Y to the nearest integers towards infinity, and D i jmax is an extreme cumulative damage value for the elements group with the i-th design variable.If the new point , the following correction has to be realized: The presented approach assumes that the design variables are arranged in ascending order.Based on the described approaches, we developed our computational program DISCRET_OPT_FAT (in MATLAB) with the fundamental algorithm and methodology presented in Fig. 3.

Optimizing discrete thicknesses in the case of the thin-shell finite element model
Let us design I cross section of the curved beam, as in Fig. 4, which is excited by random forces Fx, Fy, Fz (see Figs. 5-8).The design variables will be the thicknesses X1, X2, X3 (see Fig. 3).The computational parameters will be considered by the following: Young's modulus We then apply the new discrete optimizing approach and the Gauss-Seidel optimizing method for a comparison study.The optimum will be found using the parameters presented in Tab. 1.The chosen results of the optimizing process and a short comparison study are presented in Tab. 2 and Figs.9-12.