The “ Projection-by-Projection ” ( PbP ) criterion for multiaxial random fatigue loadings : Guidelines to practical implementation

This work is motivated by the increasing interest towards the application of the “Projection-by-Projection” (PbP) spectral method in finite element (FE) analysis of components under multiaxial random loadings. To help users and engineers in developing their software routines, this paper presents a set of numerical case studies to be used as a guideline to implement the PbP method. The sequence of analysis steps in the method are first summarized and explained. A first numerical example is then illustrated, in which various types of biaxial random stress are applied to three materials with different tension/torsion fatigue properties. Results of each analysis step are displayed explicitly to allow a plain understanding of how the PbP method works. The examples are chosen with the purpose to show the capability of the method to take into account the effect of correlation degree among stress components, and the relationship between material and multiaxial stress in relation to the tension/torsion fatigue properties. A case study is finally discussed, in which the method is applied to a FE structural durability analysis of a simple structure subjected to random excitations. The example describes the flowchart and the program by which to implement the method through Ansys APDL software. This final example illustrates how the PbP method is an efficient tool to analyze multiaxial random stresses in complex structures.


INTRODUCTION
t may happen that a new multiaxial fatigue criterion, first published as an article in a scientific journal, gains an increasing interest also in the technical community, among engineers or developers of fatigue-based software.It may also happen, though, that the scientific article (clearly more focused in explaining the theory than providing detailed results or case studies) does not give enough practical information to allow users or engineers to implement correctly the criterion on their I own.Users, in implementing the criterion, may feel not sure to have correctly understood the basic principles of the method and have doubts that their own programs are really free from errors and/or misinterpretations, and that they return a output truly correct.The possibility to have access to detailed results for a number of simple case studies, to be used as a reference, may thus help to solve possible doubts.The "Projection-by-Projection" (PbP) method is fatigue criterion for multiaxial loadings.It has first been proposed about ten years ago [1,2] as a conventional time-domain criterion, in which stress time-histories are processed directly.It was then reformulated in the frequency-domain to address multiaxial random loadings [3,4].Indeed, like for other multiaxial criteria, a time-domain definition may become computationally time-consuming (and thus impractical), for example when considering medium-to-large finite element models, in which it is required to process long digitalized multiaxial timehistories in hundreds or thousands of nodes [5].This is perhaps the main limitation (along with others theoretical issues here not mentioned) that motivated the researchers' effort of reformulating multiaxial criteria (PbP included) from time-to frequency-domain.Unlike their time-domain counterparts, frequency-domain multiaxial criteria (also called multiaxial spectral methods) are indeed able to reduce the overall computational time drastically, while keeping high levels of accuracy.In multiaxial spectral methods, the fatigue damage and life are estimated directly from Power Spectral Density (PSD) data (i.e.spectra and crossspectra), which characterize a multiaxial random stress in the frequency-domain [5].This approach proves to be computationally efficient especially if it can exploit results of FE structural dynamic analyses.It has furthermore been demonstrated that the PbP method, unlike other frequency-domain criteria, is sensitive to the local state of stress in a structure as it correctly takes into account the degree of correlation between stress components.At the same time, the method can also account for different material behaviors depending on different types of multiaxial stress [1][2][3][4].This is of great relevance in engineering applications, in which several materials need to be scrutinized and compared for the same component geometry, or it is required to perform a geometric optimization to make the best use of a given material.Though the PbP spectral method has been presented in a number of scientific articles published so far [3,4], it is our belief that this work -entirely devoted to its practical application -would be helpful to anyone who wishes to implement the method on his own.After summarizing some theoretical formulas to characterize deviatoric and hydrostatic stress components in frequency-domain, the work discusses all the steps necessary to implement the PbP method.A first numerical case study is considered, in which various biaxial stress states are applied to three different materials.Each analysis step is commented on in detail to give an exhaustive description of the procedure.Finally, the case of a simple 2D structure subjected to random excitations is considered to show how the PbP method can be easily embedded into a FE durability analysis.

FREQUENCY-DOMAIN DESCRIPTION OF A MULTIAXIAL RANDOM STRESS
efore discussing the PbP criterion and the numerical examples, it is useful to summarize the main quantities used in the following (a nomenclature is also given at the end of the article).The next paragraphs will assume that the reader is somehow familiar with the basic concepts and terminology of the frequency-domain approach to random processes.Otherwise, a short review is provided in Appendix A. Equations will be developed for a biaxial stress state; extension to a three-dimensional stress state is straightforward.

Description of multiaxial stress state in physical space
Assume that σ(t) represents the time-varying stress tensor in a point of a mechanical component (t is the time variable).If the point is located on the surface, where fatigue cracks usually nucleate, the tensor σ(t) will have three non-redundant stress components σx(t), σy(t), τxy(t) (i.e.biaxial or plane stress), where σ and τ are normal and shear stress, respectively.Each stress σx(t), σy(t), τxy(t) is assumed to be a zero-mean and covariant stationary random stress (the term "stationary" means, roughly speaking, that the random stress has statistical properties almost constant over time).The stress components in σ(t) are usually reordered into a stress vector x(t)=(σx(t), σy(t), τxy(t)), also written as x(t)=(x1(t), x2(t), x3(t)).This vector is characterized in the frequency-domain by the two-sided PSD matrix:

S
(1) B in which diagonal terms are auto-spectra and out-of-diagonal terms cross-spectra; f is the frequency (Hz).Matrix S(f) is Hermitian (see Appendix A).
The PSD matrix allows the covariance matrix C (symmetric) to be computed: The i-th diagonal term Cii=Var(xi(t)) is the variance of stress xi(t), the ij-th out-of-diagonal term is the covariance C ij =Cov(x i (t), x j (t)) between x i (t) and x j (t).Normalizing the covariance C ij allows a correlation coefficient to be defined as . The correlation coefficient represents, for multiaxial random loadings, the statistical equivalent of the phase angle for sinusoidal loadings.It discriminates between two limiting cases: rij=1 for perfectly correlated processes (i.e.proportional stresses), rij=0 for uncorrelated processes (i.e.not-proportional stresses).

Description of deviatoric and hydrostatic stress
The PbP method is an invariant-based multiaxial criterion, so a frequency-domain description of the deviatoric and hydrostatic stresses is needed.The PbP criterion works with the amplitude a 2, J of the square root of the second invariant J2 of the deviatoric stress tensor σ'(t), where the decomposition σ(t)=σ'(t)+σH(t)I into deviatoric and hydrostatic tensors is used.The definition relates the second invariant to the stress vector s(t)=(s 1 (t), s 2 (t), s 3 (t)) in the threedimensional deviatoric space.The following transformation rules apply [6]: Expressions (3) allow both the hydrostatic and deviatoric stress components to be directly computed from the stress vector x(t) in the physical space.When the stress tensor σ(t) changes over time, the tip of vector s(t) describes a curve (loading path ) in the deviatoric space.In a time-domain approach, conventional definitions (e.g.Longest Chord, Minimum Circumscribed Circle, etc. [6]) are used to identify the amplitude a 2, J of the loading path.In a frequency-domain approach, this time-domain procedure is replaced by a spectral characterization of the stress quantities in Eq. ( 3).As a first step, the correlation matrix of s(t) is defined as: where R() is the correlation matrix of x(t), with  a time lag (see Appendix A).The (symmetric) covariance matrix of s(t) results in the special case =0: where C=R() is the covariance matrix of x(t), see Eq. ( 2).
Similarly to x(t), the PSD matrix S'(f) of vector s(t) is the Fourier transform of the correlation matrix R'() in Eq. ( 4): The previous result yields by the fact that the Fourier transform     is a linear operator and A is a matrix of constants.By following the same procedure, it is straightforward to find the PSD expression of the hydrostatic stress σH(t): Its zero-order spectral moment (i.e.area of SH( f )) gives the variance: Expressions ( 6)-( 7) characterize the deviatoric and hydrostatic stress in the frequency-domain completely.

ANALYSIS STEPS OF THE PBP CRITERION
his Section summarizes the main analysis steps to be followed when implementing the PbP criterion (see Fig. 1).If the PbP method is applied to the output of a FE analysis, the steps have to be repeated for each nodal result in the model.The input data required by the analysis are:  PSD matrix S(f) of the biaxial stress (σ x (t), σ y (t), τ xy (t)), along with the covariance matrix C in Eq. ( 2) calculated from S(f).
The stress may refer to a physical point in a structure, or to a node in a FE model (in which case matrix S(f) is output directly by a FE analysis).If S(f) is not known, it may be estimated from measured time-histories;  parameters of tension and torsion S-N curves: strength amplitudes σA, τA at NA=210 6 cycles and inverse slopes kσ, kτ.
The S-N lines may represent design curves at prescribed survival probability (97.7%).
Figure 1: Analysis steps and quantities involved by the PbP spectral method in frequency-domain.
Step 1 -From stress vector to deviatoric/hydrostatic stress The first step is to characterize the stress vector s(t) by its PSD matrix S'(f) in Eq. ( 6), and by its covariance matrix C', see Eq. ( 5) or equivalently compute C' from S'(f).Then, characterize the hydrostatic stress σ H (t) by its power spectrum S H (f) in Eq. ( 7), and by its variance VH, see Eq. (8).

INPUT
stress PSD covariance matrix

S(f), C
Material properties σ A , τ A (at N A cycles) k σ , k τ

ANALYSIS STEPS OUTPUT
Deviatoric/Hydrostatic Step 2 -From the original coordinate system to the (rotated) principal coordinate system In general, matrix C' is not diagonal, i.e. some correlation exists between the stress components of s(t).If this happens, it is necessary to compute a new covariance matrix C'p in which all cross-correlations (out-of-diagonal elements) are zero.This outcome yields by solving the following eigenvalue/eigenvector problem for C': The eigenvalues are in the main diagonal of C'p, the eigenvectors are the columns of U.
In a time-domain perspective, the eigenvalue/eigenvector problem (9) can be viewed as the projection of the loading path  along the axes of a new (rotated) "principal coordinate system".Fig. 2 shows an example for a tension/torsion loading.The non-null deviatoric stress components s 1 (t), s 3 (t) give the twodimensional loading path .Once projected in the principal system this path gives rise to two stress projections Ωp,1(t), Ω p,3 (t).Both vectors s(t) and Ω(t) then share the same loading path in the deviatoric space.
Figure 2: Example of random loading path  in the two-dimensional deviatoric space, resulting from a tension/torsion loading.The rotated "principal coordinate system" is located by axes (s 1,0 , s 3,0 ) and it is used to obtain the stress projections Ω p,1 (t), Ω p,3 (t).
In a frequency-domain approach, the direct projection  is yet not necessary, as it is replaced by a "rotation" of the PSD matrix from the "old" to the "new" coordinate system.Indeed, vector Ω(t) is characterized in the frequency-domain by the PSD matrix S'p(f), which turns out by "projecting" the spectral matrix S'(f) in the new rotated principal system: 10 00 15 00 20 00 25 00 Matrix S' p (f) is diagonal, by definition, and it collects the auto-spectra S' ii (f) of stress projection Ω p,i (t), i=1, 2, 3. Thanks to the fact that, in the principal system, the stress projections Ω p,i (t) are completely uncorrelated, their fatigue damage can be computed separately.Knowing the power spectra S'pii(f), the damage can be estimated through spectral methods for uniaxial loadings (see Step 4).To this end, each spectrum needs to be characterized by several parameters (see Appendix A): spectral moments (λ 0i , λ 1i , λ 2i , λ 4i ), bandwidth parameters (α 1i , α 2i ), frequency of up-crossing and peaks (ν 0,i , ν p,i ). Step

-Reference S-N line in the Modified Wöhler Diagram
The fatigue damage for each projected stress Ωp,i(t) is calculated on a "reference S-N curve" in a Modified Wöhler Diagram (MWD) [7].This diagram relates the number of cycles to failure, N, to the amplitude of the square root of second invariant of stress deviator, (from now on, this simplified notation will be used to avoid the square root symbol).

  
and k τ >k σ .Other arrangements of S-N lines are yet possible depending on the combination of fatigue properties characterizing a specific material.The list in Tab. 1 shows that materials are indeed characterized by S-N properties over a wide range [8].

Material
Ref.The reference S-N line in the MWD is located by the stress ratio [3,4]:

Type of loading σ
which is a function of the mean value σH,m and variance VH of the hydrostatic stress σH(t), and of the square root of the sum of variances C p,ii of stress projections Ω p,i (t).Parameter ρ ref only depends on the multiaxial stress, not on material.In Eq. ( 11), the numerator approximates the maximum of hydrostatic stress, whereas the denominator estimates the equivalent amplitude of cycles counted in each projection.If all stress components have a zero mean value, it is σH,m=0.In case of constant amplitude multiaxial loading (sinusoidal), the expression (11) returns the definition of ρ ref given in [1,2] for the time-domain criterion.
Parameter ρ ref quantifies the relative contribution of hydrostatic to deviatoric stress components in a multiaxial stress.Two limiting cases exist for uniaxial loading: a stress ratio ρref=1 for tension or bending (only normal stress), ρref=0 for torsion (only shear stress).A purely hydrostatic state of stress would have ρ ref .On either two limiting cases (i.e.ρ ref =0 or 1), the reference S-N line coincides with the line of the corresponding uniaxial loading (tension or torsion).In any other case in which the loading is multiaxial (i.e. for any other value 0≤ρref≤1), the reference S-N curve would lie between those for tension and torsion, its position being established by ρ ref .
from the amplitude strengths JA,, JA,τ and inverse slopes kσ, kτ of the tension and torsion S-N curves.
Step 4 -Fatigue damage calculated for each stress projection Each stress projection Ω p,i (t), obtained in Step 2, is a uniaxial random stress.Its damage in time unit (damage/s), say d(Ω p,i (t)), can be estimated in the frequency-domain by uniaxial spectral methods.No restriction is imposed on which method to use from those available in the literature [15,16], although "wide-band methods" are recommended if the PSD of each projection Ω p,i (t) is not narrow-band.For example, fairly accurate estimations are given by the "Tovo-Benasciutti (TB) method" [15,16]: where Γ(-) is the gamma function and η TB,i is a correction factor that accounts for the spectral bandwidth of S p,i (f) and it depends on a proper weighting coefficient bapp (for their expressions, see [15,16]).In the limiting case of narrow-band PSD, η TB,i →1.Eqn.(13) has been proved to provide estimations close to those from Dirlik's method, which the reader may be more familiar to [17]: Parameters D 1,i , D 2,i , D 3,i , Q i , R i are best-fitting coefficients; their expressions can easily be found in the literature (see, for example, [15][16][17]).
J a (log) Step 5 -Estimating the total fatigue damage Once the fatigue damage d(Ω p,i (t)) for every projection has been calculated, the total damage d(Ω) caused by stress vector Ω(t) (which coincides with the damage of the multiaxial stress x(t)) can be estimated by the following expression: This expression employs a non-linear combination rule to handle damaging events in uncorrelated stress projections.This law is able to account for the phase shift between stress components.Its validity versus experimental data has been checked in previous publications [3].Eqn. ( 15) is very general.Its mathematical expression depends on the particular spectral method used for estimating d(Ωp,i(t)).If the "TB method" is considered: The Dirlik's method yields, instead, a slightly more complex equation: The

Test cases with biaxial stress
o further clarify the analysis steps sketched in Fig. 1, the PbP method is now applied to some simple case studies, in which 4 types of biaxial stress are combined with 3 types of material fatigue properties.The decision to use simple case studies is no coincidence.In fact, it permits a much clearer monitoring of results in each analysis step, which in turn makes the understanding of the PbP method much easier.On the other hand, the case studies here analyzedthough simple -do synthesize combinations of materials and stress states that are actually encountered in mechanical components subjected to multiaxial loadings.The examples, in particular, are also specifically conceived so to better emphasize some peculiarities of the PbP method.Each biaxial stress has components σ xx (t), σ yy (t), τ xy (t) that are zero-mean stationary and Gaussian, and with a band-limited rectangular PSD, see Fig. 4. Also the use of a rectangular PSD, in place of other more realistic (but irregular) spectra, contributes to make results much easier to be understood.The rectangular one-sided PSDs are centered on frequency fc=30 Hz and are characterized by a fixed maximum-to-minimum ratio fmax/fmin=15, which guarantees that PSDs are wide-band with α 1 = 0.893, α 2 = 0.771 (these parameters remain unchanged for the power spectra of deviatoric and hydrostatic stress, T as we as for stress projections).Being the frequency range fixed, each PSD is fully defined by the height h of each rectangle.Through the rectangle area, the height controls the variance Vii=hi,i•(fmax-fmin) of auto-PSDs and the covariance Cij=hi,j•(fmaxfmin) of cross-PSDs (subscript i and j stands for the stress component xx, yy and xy).It is easy to derive the relationship between the cross-PSD height and the correlation coefficient r i,j .The height h i,j is then bounded as , see Fig. 4. The cross-PSD height varies according to the different correlation degrees in the range ri,j=01.Three hypothetical material models, with different fatigue properties, are investigated:  Material A: this material has tension and torsion S-N lines with identical slope (k=kτ) and with amplitude strengths σA and τA=σA/√3, i.e. the torsion line is scaled by 3 .In the MWD, the tension and torsion S-N lines overlap each other, and they also overlap the reference line for any value of ρref.This material model has specifically been conceived for being not sensitive to the type of applied multiaxial stress (i.e. the material is not sensitive to ρ ref ).This situation reflects the hypothesis that the material obeys the Von Mises equation under fatigue loadings, and it is only sensitive to the deviatoric stress component while not sensitive to hydrostatic stress;  Material B: this material model has tension and torsion S-N lines with identical slope (k  =k τ ), but with fatigue strengths not scaled as in von Mises criterion, that is τA≠σA/√3.Unlike Material A, this material has a fatigue strength that depends on the hydrostatic stress;  Material C: this material model has arbitrary tension and torsion S-N lines, i.e. fatigue strengths as per Material B, but with different slopes for tension and torsion (k≠kτ).This is the most general situation.Like Material A, this material is sensitive to the hydrostatic stress, too.In Tab. 3, the tension S-N curve is kept fixed, while only the torsion line moves upward from case A to C. Table 3: Fatigue properties of three different material models, considered in numerical simulations.
Tab. 4 gathers the results of all the load cases examined.In particular, it collects some quantities calculated throughout the analysis: the non-zero elements of covariance matrix C' in deviatoric space and the variance of hydrostatic stress (Step 1), the non-zero elements of covariance matrix C' p in the principal coordinate system (Step 2), the parameters of the reference S-N curve in the MWD (Step 3), the damage of each stress projection (estimated by narrow-band formula and TB method) (Step 4), the total damage estimated by PbP method using the narrow-band formula or TB method (final analysis Output).
Load Case 1A and 2A refer to biaxial normal stresses applied to material A; Case 1A has "in-phase" stresses (r xx,yy =1), Case 2A "out-of-phase" stresses (rxx,yy=0).The comparison of damage values shows that Case 2A is more damaging than Case 1A.This result is due to the greater contribution of deviatoric stress components in Case 2A compared to Case 1A (see Step 2 in Tab.4), whereas it does not depend on the hydrostatic stress component, which is higher in Case 1A.Though the hydrostatic stress affects the damage by changing the S-N reference curve via ρref value, Material A has specifically been conceived for no reaction to a change in the hydrostatic stress.The comparison Case 1A vs. 2A shows that not-proportional stress lead to a greater damage.This effect of stress correlation cannot be generalized, yet.It closely depends on the material type.For example, Tab. 4 makes clear that a different material (B or C) would behave differently.For example, Case 3B and 4B highlight that damage in Material B is not correlated to the proportionality degree.
In Tab. 4, it is worth noting that Cases 2A and 3A yield the same damage value.This outcome may be explained by considering the role of deviatoric and hydrostatic stress components.On one side, either Cases have in common (see Step 2 in Tab. 4) the sum (C' p22 +C' p33 ) of variances of deviatoric stress components in the principal coordinate system (the third projection has C'p11=0).This sum represents a sort of "total variance" of stress projections Ωp,2(t), Ωp,3(t) and it enters the damage expression ( 16) in the special case in which the bandwidth parameters η TB,i , ν 0,i take on identical values for all stress projections (a case that actually occurs in all stress states of Tab. 2).It is trivial to verify that damage d(Ω) is proportional to  Predictably, the trends in Fig. 5 confirm how different materials respond differently to the same multiaxial loading.While materials B and C show a moderate-to-great sensitivity to the type of loading, material A shows no sensitivity at all.In particular, Fig. 5 makes crystalline clear that material A is unable to discriminate between different multiaxial states of stress; the estimated damage is constant and always equal to the value for torsion loading.Materials like A, however, may be very common in engineering applications.It may characterize those unfortunate situations in which only the tension S-N curve is available from experiments, whereas the torsion S-N curve is only approximated by taking a line parallel and moved downward by 3 (Von Mises rule).Although not recommended, this A-type material represents the only solution feasible for applying the PbP criterion, if no fatigue data for torsion loading are available.It has to be emphasized that such a "Von Mises approximation" is even hypothesized (implicitly) in some multiaxial criteria (e.g.[18]), which may thus lead to unrealistic estimations (see [4,8]).

FE durability analysis an L-shaped structure under random acceleration
The purpose of this second example is to illustrate how the PbP method can easily be embedded into a FE durability analysis of a structure undergoing random vibrations.The example considers an L-shaped structure made of steel, whose geometry imitates that one already studied in [18]; some dimensions have been slightly changed [19] to enhance stress concentration effect at the hole and two lateral notches.A finite element model with "shell" elements is used to discretize the structure.A mapped mesh is generated in the notch regions.The average elements size is 3.5 mm; the smallest element dimension of 0.8 mm appears in regions of mesh refinement at notch tip.The model has a total of 1529 elements and 1687 nodes.The structure is clamped at both ends, at which a random acceleration is imposed along the direction normal to the specimen plane.Input accelerations have a band-limited (rectangular) one-sided PSD, ranging from 1 to 200 Hz, with height 25π (m/s 2 ) 2 •Hz -1 (as in [18]).Input accelerations applied at the two clamped ends are fully correlated (their cross-PSD is different from zero).A frequency-domain spectrum analysis is carried out through FE simulations to determine the structure natural frequencies and the stress PSD matrix in every FE node.Stress spectra are next processed by the PbP method to determine the total fatigue damage caused, in each node, by the local multiaxial state of stress.The whole numerical analysis is performed by software Ansys with APDL language, which is also used to implement the PbP method.Appendix B provides more details on the numerical algorithm.This analysis strategy, in which all calculations are carried out in Ansys software, simplifies the overall data management and thus represents an enhancement of the two-phase procedure followed in [4], in which nodal FE results (PSD matrices) were first calculated in Ansys, then exported and post-processed in Matlab to apply the PbP criterion.The first five natural frequencies returned by modal analysis are: f1=16.1 Hz, f2=67.3Hz, f3=82.2Hz, f4=175.7 Hz, f5=178.6Hz.
In each node the state of stress turns to be biaxial.Fig. 7 displays an example of auto-and cross-PSDs in a node close to the round notch.Results refer to the "top" layer of shell elements.For fatigue damage estimation, two different types of material (Material 1 and 2) were considered for the structure.Material 1 (likewise material A in Tab. 2) is only sensitive to deviatoric stress components, whereas Material 2 (likewise material C in Tab. 2) is also sensitive to hydrostatic stress component.The material properties used in simulations are: σA=300 MPa,      Likewise, the trend in Fig. 9(b) remarks how much the fatigue damage (via PbP criterion) is influenced by the hydrostatic stress via ρ ref values.Being dependent on hydrostatic stress, Material 2 is characterized by damage values that depend upon ρref more clearly and that differ by up to an order of magnitude from those characterizing Material 1.Though simple, the numerical example discussed so far shows, in fact, the ability of the PbP method to capture the material sensitivity to the local multiaxial stress in a structure.This result is of great relevance in engineering applications, especially at early design phases, when it is desirable to compare the performance of different materials for the same component or, by contrast, it is required to identify a material whose fatigue properties are the most suitable for a given component subjected to a known vibration input.

CONCLUDING REMARKS
he writing of this article was primarily motivated by the aim to provide engineers engaged in vibration fatigue with a practical guide to apply the PbP spectral method.After a summary of the main analysis steps of the PbP criterion, the article applied the criterion to some simple numerical case studies, in which three material types are subjected to four different types of biaxial random stress (e.g.tension-tension, tension-torsion, "in-phase" or "out-of-phase").A simple rectangular PSD for the stress is considered to allow a much easier understanding and interpretation of the obtained results.The results of each analysis step were listed in detail to serve as a benchmark for users interested in implementing the PbP method by their own.Results are also used to highlight the main features characterizing the PbP criterion.For example, the examples allow the role of material sensitivity to hydrostatic and deviatoric stress to be clearly pointed out.The case studies make apparent the capability of the PbP method to take into account the correlation degree between stress components, as well as the different material behavior as a function of the various states of stress.The article concludes by an example in which the PbP method is applied to the FE-based fatigue analysis of a thin structure submitted to a random acceleration.This example shows the advantage of the PbP method in CAE-based design of structural components undergoing multiaxial fatigue loading.It emphasizes, in particular, how the method proves to be an efficient tool in the design of complex structures, e.g. for discriminating among materials with different fatigue properties, or for identifying the material with the most suitable properties.
In practical applications, where negative frequencies have no direct physical meaning, the two-sided spectrum S(f) is replaced by a one-sided spectrum limited to positive frequencies only, which is defined as G(f)=2S(f), 0<f<∞ and zero elsewhere.It is customary to describe S(f) or G(f) by the set of spectral moments [21]: Eqn. (19) shows that the variance of x(t) is the zero-order moment Var(x(t))=R(0)=λ0, which corresponds to the area of G(f).If x(t) is Gaussian, the frequency of zero up-crossings, ν0, and the frequency of peaks, ν0, are [21]: Spectral moments are also combined into bandwidth parameters, as for example [21]: where 0≤α m ≤1 and α 1 ≥α 2 .Bandwidth parameters summarize the shape of a PSD.Two limiting cases exist: a narrow-band PSD with α1→1, α2→1, a wide-band PSD with α1<1, α2<1.The quantities in Eq. ( 19)-( 21) enter the analytical expressions used by spectral methods for estimating the fatigue damage.
The previous definitions can be generalized to a biaxial stress x(t)=(σ x (t), σ y (t), τ xy (t)).By analogy with Eq. ( 18), x(t) is characterized in time-domain by a correlation matrix R()=E[x(t)•x(t+)] and in frequency-domain by a PSD matrix, which are a Fourier transform pair:

) and Rij()=E[xi(t)•xj(t+)] (i≠j) is the cross-correlation functions between xi(t) and xj(t). In general Rji()≠Rij().
The PSD matrix takes the form: Diagonal elements are auto-PSDs, out-of-diagonal elements are cross-PSDs.Each auto-PSD S hh (f) is a non-negative realvalued even function of f.Cross-PSDs are complex-valued functions of f and can be written as , which implies that Shk(f) must be zero at any frequency where either S hh (f) or S kk (f) is zero too, whereas it maximum depends on the values of the auto-PSDs.Also in the multiaxial case, a one-sided PSD matrix can be introduced: G(f)=2•S(f) for f>0, zero elsewhere.Similarly to Eq. ( 19), a matrix of spectral moments of S(f) is defined as: Diagonal elements λn,hh are the moments of auto-PSDs, out-of-diagonal elements λn,hk (h≠k) the moments of cross-PSDs: Eqn. (25) shows that the cross spectral moment λ n,hk is only function of the co-spectrum ) ( f S c hk , as the quad-spectrum ) ( f S q hk is an odd function of f and it simplifies in the integral computation (25) (this finding has somehow to be expected, since spectral moments are real numbers).
Eqns. ( 22) and (24) show that =0 yields the covariance matrix C, which also coincides with the matrix λ 0 of zero-order moments.Variance and covariance terms are: This final result also confirms how cross-PSDs (which are complex-values functions) are not involved in the covariance matrix (which is real-valued).

APPENDIX B -IMPLEMENTATION OF PBP METHOD IN ANSYS APDL
ig. 10 displays the flowchart of the PbP method, along with the scripts used to implement the method in Ansys software.The relationship with the previous approach in Matlab is also shown for comparison.In both approaches, the analysis is developed in two separate phases that must be executed in sequence.The first phase is managed in Ansys by a main program that generates the finite element model (mesh and boundary conditions) and computes the stress PSDs in each node ("DO/ENDO" loop) through the procedure of random vibration analysis.Throughout this phase, the stress PSD for each node is stored into a text file (function "psdcovdata2text.mac").All text files are gathered in a subfolder for subsequent processing.Similarly, also the finite element model is saved into a *.cdb file.The list of nodes and elements need to be exported and stored (function "FEnode_elem.mac")only if the next analysis phase is performed in Matlab.The second phase, carried out either in Ansys or in Matlab, computes the damage in each node according to the PbP method.In Ansys, this phase is managed by a main program through a "DO/ENDO" loop, in which the stress PSD data for each node are first retrieved from text files (function "text2psdcovdata.mac")and then used as input for subsequent analysis.Also the S-N material properties are specified at the beginning of the analysis.At the end of second phase, results are displayed as a contour plot (damage map), or written into a text file (function "PbPresults2text.mac").The first phase must be run one time only.Once stress PSDs have been calculated and stored, the second analysis phase is performed to compute the fatigue damage.This phase is repeated as many times as are the combinations of S-N properties that need to be scrutinized.In terms of computation time, the first phase may be slightly longer than the second one (especially with low RAM memory).Just to provide an order-of-magnitude estimate, for the model in Fig. 6, determining the nodal stress PSDs takes about 4 minutes, while applying the PbP method takes about 3 minutes on a 64-bit workstation (CPU 3.80 GHz, 32 GB RAM).Compared to the Ansys/Matlab mixed procedure, the new approach entirely based on Ansys has the great advantage to exploit the graphical capabilities of the commercial software in displaying the damage contour maps, especially in 3D complex finite element models.Another little advantage, though trivial, is that only one type of software is needed to carry out the whole calculation.On the other hand, the main disadvantage of using Ansys as the only computation tool is due to a greater complexity and less flexibility of APDL language in performing the fatigue damage calculations required in the PbP method.Although the APDL language -as Matlab -can execute a lot of algebraic and trigonometric functions, and it can even manage "DO/ENDO" and "DOWHILE/ENDO" loops, it does not have functions to compute the eigenvalues/eigenvalues of a square matrix or the gamma function Γ(-) directly.Therefore, an APDL macro (called "eigen.mac")has been written to deal with eigenvalues/eigenvalues computation.This macro took the cue from the software library EISPACK, written in Fortran [22,23].In particular, the functions involved are: TRED2 for tridiagonalization of a symmetric matrix [24], TQL2 for the QR eigenvalue algorithm [25].Instead, the gamma function Γ(-) was computed by means of Stirling's approximation formula, see for example [26].To conclude, it is worth emphasizing that the flowchart here described for Ansys software is very general and it can easily be adapted to any other finite element code that perhaps the reader knows better.

Fig. 3
Fig. 3 depicts the relationship between the MWD and the Wöhler diagram.The figure shows, on the right, the reference S-N line along with the tension and torsion S-N lines in the MWD, and it also clarifies how the tension/torsion lines in MWD have to be sketched from the corresponding lines in the Wöhler diagram on the left.From the first and third expression in Eq. (3), in the MWD the reference strength amplitudes at NA cycles are 3 / , A A J

Figure 3 :
Figure 3: Relationship between S-N lines in Wöhler diagram (left) and Modified Wöhler Diagram (right).The reference S-N line for a general multiaxial stress (ρ=ρref) is also shown.Symbol Ja stands for a 2, J .In a log-log diagram, the reference S-N line is expressed by the equation

Figure 4 :
Figure 4: Example of one-sided auto-and cross-PSD for a tension-torsion loading with non-zero correlation rxx,xy (stress case 3).

.
On the other side, damage (16) is computed on the same S-N curve for both Case 2A and 3A, despite they have a different variance V H for the hydrostatic stress and, hence, a different ρ ref value (see Tab. 4).However, as already said, for Material A the hydrostatic stress has no effect on damage.Cases 3A, 3B and 3C emphasize the effect of different materials under the same tension-torsion loading (stress case 3).It is not surprising that a change of material determines a change in the estimated damage, too.This result is governed by the reference S-N line, which takes on different positions for each material type, as confirmed by the dissimilar values of parameters J A,ref , k ref (see Step 3 in Tab. 4).For example, the damage increases from material A to C, in relation to the different positions assumed by the reference S-N line.

Figure 5 :
Figure 5: Effect of hydrostatic stress on fatigue damage for different materials undergoing a tension-torsion loading with C s1 +C s3 = 1.The graph in Fig. 5 allows the role of material properties to be further pointed out.The figure depicts the trend of the fatigue damage d TB (Ω) in the case of a tension-torsion loading, as a function of ρ ref in the range 01.Parameter ρ ref synthesizes different types of tension-torsion multiaxial loading, bounded by the limiting uniaxial case of pure torsion (ρref=0) to pure tension (ρ ref =1).The PSD of normal and shearing stress are suitably normalized so to keep constant the variance sum C s1 +C s3 =1, as well as the corresponding sum const C ii p    ,

Figure 6 :
Figure 6: Geometry and finite element model of the L-shaped beam.Thickness is 0.5 mm.

Fig. 8
displays the damage maps for either material.The damage is estimated by PbP-TB method in Eq. (16) and refers to the "top" layer of shell elements.Left Figs.(a) show the damage D in linear scale to better highlight the locations of the most damage points; right Figs.(b) show the logarithmic values log 10 (D) to allow the differences in the overall damage distribution to be appreciated best.The comparison of figures makes evident how the damage distribution does change with material type.In particular, the location of the most damaged point in the structure shifts in accordance to the degree by which material is sensitive to the hydrostatic stress value.

Figure 7 :
Figure 7: Auto-and cross-PSD (one-sided) in a node close to the round notch.

Figure 8 :
Figure 8: Comparison of damage maps for Material 1 vs. Material 2. Figs.(b) plots of log 10 (D) from Figs. (a).Indeed, in a FE model the stress ratio ρref may vary over a range of values, which depend on the variation of the local multiaxial state of stress among nodes.Fig. 9(a) illustrates the distribution of ρ ref values obtained in numerical simulations.This distribution holds for both materials, being ρ ref only a function of the local stress, not of material.The figure also marks the values at which the multiaxial stress is close to uniaxial tension or torsion, and the regions A, B, C, D identified in Fig. 8(a) (top).The values of ρ ref then measure if the stress in a point is more tension-like or torsion-like.The capability of the PbP criterion to capture the material sensitivity to the local stress state is well summarized by the graph in Fig. 9(b), which plots the values of the damage ratio DMat2/DMat1 vs. ρref for each node of the FE model.If the PbP were not sensitive to the material which the structure is made of, the FE analysis would return the same damage value in every node regardless of the material type, and in turn all the dots in Fig. 9(b) would lie on the dashed horizontal line.

Figure 9 :
Figure 9: (a) Distribution of ρ ref values in FE nodes; (b) damage ratio D Mat2 /D Mat1 vs. ρ ref values (each point refers to one nodal result).
t) be a zero-mean uniaxial random stress.It is characterized, in time-domain, by the autocorrelation function R()=E[x(t)•x(t+)] ( is a time lag) and, in frequency-domain, by a two-sided Power Spectral Density (PSD) function S(f), ∞<f<∞.Both functions constitute a Fourier transform pair (Wiener-Khintchine relations)[20]: i is the imaginary unit.The real part -spectrum") is a real-valued even function of f, -spectrum") is a real-valued odd function of f.The cross-PSD can be positive or negative.Matrix S(f) is Hermitian because means complex conjugate.Finally, the values of the cross-PSD are also bounded by the Schwartz'

Figure 10 :
Figure 10: Flowchart of the APDL script to implement the PbP in Ansys software.The relationship with the previous Matlab-based procedure (sketched on the right) is shown.
covariance matrix of x(t) and s(t) C'p covariance matrix in the principal coordinate system dTB(Ωp,i(t)) damage of stress projection Ωp,i(t) by TB methodd(Ω) total damage for stress vector Ω(t) E[-] expected value G(f) one-sided PSD matrix of x(t) k σ , k τ inverse slope of tension and torsion S-N curve square root of second invariant of stress deviator J A, , J A,τ , k  , k τ amplitude strengths and inverse slopes of the tension and torsion S-N curves in MWD J A,ref , k ref amplitude strength and inverse slope of the reference S-N curve in MWD NA reference number of cycles r ij correlation coefficient between x i (t) and x j (t)

Table 1 :
Fatigue parameters for several plain materials or notched specimen).

Table 2 :
(1,2,3,4) the stress cases(1,2,3,4)in Tab. 2 with the materials models (A, B, C) in Tab. 3 results in the whole set of load cases (1A, 2A, 3B, …) in Tab. 4, which are examined in the numerical examples.The stress cases in Tab. 2 are: All through the examples, the normal stresses have variance V xx =V yy =3 and the shearing stress V xy =1, so that all components of the stress deviator share the same unitary variance.Variance and correlation degree characterizing the biaxial stress states considered in numerical simulations.

Table 4 :
Summary of step-by-step results (only non-zero values) of numerical case studies.In Step 4 and Output, the damage values are