An efficient procedure to speed up critical plane search in multiaxial fatigue: Application to the Carpinteri-Spagnoli spectral criterion

A more efficient procedure is proposed to speed up the Carpinteri-Spagnoli (CS) algorithm in numerical computations. The goal is accomplished by deriving the exact solution for the spectral moments and expected maximum peak of normal/shear stress in any rotated plane orientation. The algorithm then avoid the use of “for/end” loops to identify the five rotations that locate the critical plane in CS method. The procedure is especially advantageous if applied to three-dimensional finite element analysis, in which the stress spectra in thousands of nodes need to be processed iteratively. The procedure is based on theoretical results that have, however, a more general validity, being applicable to any multiaxial criterion that makes use of angular rotations to identify the critical plane.


Introduction
Multiaxial spectral methods are a special class of fatigue criteria which characterize a multiaxial random stress in the frequency-domain through its Power Spectral Density (PSD) matrix [1]. This article, in particular, deals with the Carpinteri-Spagnoli (CS) spectral method, in which a sequence of five rotations ( Fig. 1) is used to identify the critical plane onto which an equivalent stress is finally defined [2][3][4].
If, to implement the CS method, a routine is written by following exactly the algorithm as described in several literature articles [2][3][4], it seems that the outcome is a numerical code that is not very efficient and thus slow. This occurs for two main reasons. First, searching for the critical plane requires that two or more angles be scanned in the threedimensional space. In numerical routines, where angle ranges are digitalized in discrete values, this task can be accomplished only through a sequence of nested "for/end" loops, which make the total number of iterations to increase considerably. Secondly, at each iteration in which a given rotated plane is analyzed, the CS algorithm asks to compute a certain number of quantities (e.g. cross spectra and related spectral moments in the rotated PSD matrix) that, in fact, are not really necessary to identify the critical plane.
The computational time may grow from a few minutes up to hours or days (thus far beyond any practical limit) when each simulation run needs to be iterated hundreds or 1 Corresponding author: denis.benasciutti@unife.it thousands of times to analyze the nodal results in a three-dimensional finite element model. This is a severe limitation if the method is to be routinely applied in design.
Having this situation in mind, this works proposes a novel procedure for implementing CS criterion, which permits the computation time to be shortened significantly. The procedure takes advantage of the exact solution of spectral moments and expected maximum peak of normal/shearing stress in any rotated reference frame. This solution is the basis of a numerical routine that, being free from any "for/end" loop, is definitively much faster than the standard CS algorithm.
It has to be remarked that the theoretical approach, here devised for the CS criterion, has a validity yet more general, being it applicable to any other multiaxial spectral method that makes use of angular rotations to identify the critical plane or the equivalent stress.
A multiaxial stress state is conveniently represented by the random stress vector The zero-order moment matrix (m=0) is the covariance matrix 0,XYZ = XYZ of In what follows, also the second-order moment matrix 2,XYZ will be used. Subscripts specify the reference frame in which the PSD matrix and its spectral moments are computed; other reference frames will indeed be introduced in the following sections.

The Carpinteri-Spagnoli spectral method
The stress vector XYZ ( ) and its PSD matrix XYZ ( ) are defined in the fixed reference frame XYZ with origin in P, see Fig. 1(A). A rotated frame X ′ Y ′ Z ′ with same origin P is defined through three Euler angles ( , , ). In the rotated frame, the stress vector where ( , , ) = ψ θ ϕ is the product of three rotation matrices [2][3][4]: Throughout the text, trigonometric functions are abbreviated x = cos , x = sin . For X ′ Y ′ Z ′ ( ), the matrix of m-th order spectral moments m,X ′ Y ′ Z ′ is defined as in Eq. (3). The first step in the CS method is to scan the angles and (in intervals 0 ≤ ≤ 2 and 0 ≤ ≤ ) in order to find that particular direction ′ (defined by * and * ) which maximizes the expected maximum peak of normal stress z ′ ( ) in time [2][3][4]: where 0,3 ′ 3 ′ is the variance of z ′ (t) and 0,3 ′ 3 ′ = 0 + = √ 2,3 ′ 3 ′ 0,3 ′ 3 ′ ⁄ is the number of mean upcrossings in time length T (note that the actual value of T is irrelevant in determining the solution angles). Multiple solutions ( * , * ) may exist due to periodicity of trigonometric functions in matrix ( , , ).
Starting from X ′ Y ′ Z ′ , another frame uvw is finally identified through angles and , see Fig. 1C. Frame uvw is attached to the critical plane at point P. The stress vector in uvw T and has PSD where ′ ′ ′ * ( ) has been determined at previous steps, whereas the rotation matrix ̃( , ) = 6 represents a sequence of two rotations through matrices δ and γ , which are functions of and (the identity matrix 6 corresponds to no rotation) [2][3][4]: The off-angle δ, which denotes a clockwise rotation about the 2-axis (see Fig. 1 that only depends on normal and shear stress fatigue limits, af−1 and af−1 under fully-reversed loading [2][3][4]. The last angle to be identified is , which represents a counterclockwise rotation about the w-axis, Fig. 1 After this last step, the critical plane is fully localized by the angles ( * , * , * , * , * ). The PSD matrix of X ′′ Y ′′ Z ′′ ( ) follows directly from the PSD matrix in the initial reference frame X ′′ Y ′′ Z ′′ * ( ) =̃( * , * ) ( * , * , * ) XYZ * ( ) ( * , * , * ) T ̃( * , * ) T (the asterisk reminds that the PSD matrix no longer depends on rotation angles).
The available references [2][3][4] do not report any guideline on how to translate the CS method into a numerical algorithm. If one follows exactly the procedure as detailed above, it seems that -for any set of angles to be scanned -it is required to compute the rotated PSD matrix and its corresponding spectral moments, from which to extract only those moments ( 0,3 ′ 3 ′ , 2,3 ′ 3 ′ , 0,4 ′ 4 ′) that are necessary to determine the maximum variance or peak characterizing that particular set of angles selected previously.
In numerical algorithms, in which angle intervals are discretized, this procedure must be iterated over all the discrete values at which every angle interval is subdivided. This task is accomplished through "for/end" loops in which, at each iteration, an angle value is assigned and the resulting maximum computed. Once all digitalized angles have been scanned and the corresponding maximum stored, the overall largest maximum variance or peak can be identified, along with its corresponding angle.

Fig. 2. Numerical code of CS method: comparison between standard and new algorithm
The numerical algorithm is sketched in Fig. 2(left). The number of iterations (i.e. number of planes to be analyzed) is tot = ϕ • θ + ψ + , where is the number of points for each digitalized angle interval, which is inversely proportional to the discretization step. A narrow step improves the resolution in locating the critical plane, but also increases the number of iterations, making the algorithm computationally expensive.

Proposed algorithm
The computational efficiency of the previous algorithm can be made much faster if "for/end" loops are replaced by the analytical expressions of those spectral moments really involved in critical plane search. By applying Eq.
The subscript specifies that has a unit value in position 3 and zeros elsewhere. Equation (12) must be solved twice to determine λ 0,3 ′ 3 ′ and λ 2,3 ′ 3 ′ , which are both required in Eq. (6). This task would be rather tedious if carried out by hand calculations, but it becomes much simple with the help of a numerical computation tool (as Matlab Symbolic). A short portion of λ 0,3 ′ 3 ′ is (the complete expression is too long to be reported here):   where ii , ij are trigonometric coefficients that only depend on angles ( , , ). Also the analytical expression of 0,3 ′ 3 ′ = √ 2,3 ′ 3 ′ 0,3 ′ 3 ′ ⁄ can easily be computed from those of λ 0,3 ′ 3 ′ and λ 2,3 ′ 3 ′ . This computational artifice permits the analytical expression of [max ′ ( )] = 1 ( , , ) in Eq. (6) to be obtained as a function of the angles ( , , ) (note that the algorithm in Fig. 2(right) adopts a different strategy based on the element-byelement multiplication operator .* in Matlab).
When Eq. (6) is being solved to locate ′ , only two values * , * need to be found ( is taken zero being it irrelevant). In Matlab, this task can be achieved quite easily. First define the vectors that discretize the angular intervals 0 ≤ ≤ 2 , 0 ≤ ≤ and use them to create a two-dimensional array (command meshgrid). Then, compute 1 ( , , ) at any element in the array. Finally, look for the largest maximum in the array (or periodic maxima) by issuing the command find. This will return the solutions * , * .
Once the first two angles are found, the next step is to apply Eq. (12) with a different auxiliary vector (unit value in position 4) to determine the spectral moment in Eq. (7) directly as 0,4 ′ 4 ′ = ( * , * , ) m,XYZ ( * , * , ) T T . A quite long analytical expression is again obtained, which is very much like Eq. (13). This expression for 0,4 ′ 4 ′ = 2 ( ) is only a function of . It is then no effort to find the value * (with proper periodicity) at which 0,4 ′ 4 ′ reaches its maximum.
The previous line of reasoning is pursued to find * , * defining uvw . Equation Note that the algorithm, sketched in Fig. 2(right), is free from "for/ends" loops.

Numerical case studies
A first numerical example deals with uncorrelated random bending/torsion with narrowband power spectrum centered at 30 Hz. The covariance matrix C has non-zero elements 0,11 = 20, 0,22 = 40, 0,33 = 60, 0,44 = 200 (MPa 2 units) and zeroes elsewhere. Material has af−1 = af−1 √3 ⁄ . Five angular steps (∆=1°, 5°, 10°, 20°, 30°) have been scrutinized with the purpose of comparing the computation time needed by both algorithms and to investigate how much sensitive is the equivalent spectrum eq ( ) to the resolution Δ used to locate the critical plane. The larger resolutions (∆=10°, 20°, 30°) were included only to make the analysis the most complete, although they are actually too coarse for engineering applications.
For every Δ, either methods (standard CS and improved) have been applied to the input PSD matrix ( ) in order to determine the rotations angles ( * , * , * , * , * ) of the critical plane and the corresponding spectrum ( ). Each simulation clocked the computation time (elapsed time) needed by each algorithm to complete the critical plane search. This time is plotted in Fig. 3(A) versus the number of planes tot scanned by the standard CS method ( tot is roughly proportional to 1 ∆ 2 ⁄ ). In graphs, the lines connecting adjacent markers are only used to best emphasize the trend and do not represent actual values from simulations. For the standard CS algorithm, the computation time spans from a few seconds to about 300 seconds. For angular steps from ∆=30° down to 5°, the time only increases of about two or three times. Instead, it sharply lengthens of about an order of magnitude (from 13 up to 300 seconds) for a change from 5° down to 1° (this decrease of one fifth makes the number of planes to grow up of 23 times -from 2847 to 66063).
For the new CS algorithm, instead, the computation time remains always very small (at most, few tenths of a second) and it is also much less sensitive to the change in angular step (the green dashed line in Fig. 3(A) is indeed not really horizontal). What is quite apparent is that the new algorithm is definitively much faster than the standard one, with a maximum time saving of about 1400 times for ∆=1°.
On the other hand, it seems not necessary to attain such a small angular resolution. At least for these numerical examples, the variance of the equivalent stress seems indeed not so much sensitive to the values of ∆. Fig. 3(B) displays the trend of eq = [ eq ( )] versus ∆ (the values are normalized to the largest maximum eq,max reached for 1°). Other types of trend (not monotonic) may also be observed when processing for input PSDs. The variance is monitored as it determines the expected fatigue damage. In Fig. 3(B), only a small difference (below 10%) is observed over the angular gaps examined. The change in variance versus ∆ is caused by the deviation of the critical plane orientation with respect to the (presumably correct) position attained at very small Δ.
To further emphasize the advantage of the new procedure over the standard one, a finite element (FE) analysis of an L-shaped structure is considered. This example is of particular interest for engineering applications, being it representative of a typical design case study in which the CS algorithm needs to be iterated over hundreds or thousands FE nodes.
The FE model here analyzed consists of "shell" elements arranged in a free mesh (Fig.  4(A)), with an average element size 6.5 mm that gives a total of 394 elements and 469 nodes. The structure is subjected to band-limited accelerations applied at the two ends. The calculated stress PSDs at each node are discretized into 1247 frequency points in the range 0≤f≤200 Hz. For other details see [8].   Fig. 2 (this approach aims to mimic a FE model with very few nodes); results on the right side ( RUNS =469) come from processing the nodal stress PSDs from FE analysis (in both cases, the same frequency discretization is used).
A linear trend in log-log scale is obtained. For the standard CS algorithm, the elapsed time can be approximated as CS ≅ ( f RUNS ) Δ 2 ⁄ (for Δ<5°), where f is the number of frequency points ( f =1200 in this example). Constant = 1.445 is calibrated on the PC under use ( differs for the new algorithm). The previous formula, suitably calibrated for two small values of RUNS , would allow a quick estimate of the computation time needed for FE models with a much larger number of nodes.
The decisive advantage of the new CS algorithm is obvious from Fig. 4(B). While the new algorithm only requires 2.2 seconds for Δ=5°, the standard one demands 3790 seconds (10 hours) -that is, 1700 times more -which is a surprisingly huge time for such a small FE model. If this time were extrapolated to the number of nodes in medium-to-large-size FE model, it would go beyond any limit acceptable in practical applications.

Conclusions
Based on our understanding of the articles accessible in the literature, the critical plane CS spectral method has initially been turned into a numerical code. Following exactly the procedure as it is explained in those articles (no specific guidelines on the numerical aspect are available), the numerical code resulted to be not very efficient and rather timeconsuming. It indeed seems unavoidable for the code to use several "for/end" loop to scan, in the three-dimensional space, all the planes before the critical plane is correctly located through five rotation angles. The number of planes increases by far if the step of discretized angular values is made smaller. This numerical inefficiency could be a severe limitation if the CS criterion needs to be applied to the stress PSD in each node of a 3D FE model. This situation motivated the attempt to develop a new CS algorithm able to make the computation time much shorter. This goal is achieved by considering in the algorithm only those spectral moments and parameters that are really processed by the CS method in the critical plane search, and by deriving their analytical expressions in closed-form as a function of rotation angles and of spectral moments in the initial (fixed) frame. This strategy then avoids that unnecessary data be calculated (e.g. PSD matrix and its moments in any rotated plane to be scanned) and permits "for/end" loops to be eliminated from the numerical code.
Several numerical examples confirmed that the new algorithm can shorten the computation time as much as 1400 times for an angular step 1°. This figure of merit gives a decisive advantage of the new algorithm if the CS method is applied to the output of threedimensional FE models, in which the algorithm need to be iterated to process the stress PSDs of hundreds or thousands of nodes.
Last, but not least, it is worth noting that the procedure described in this paper has a far more general validity, its theoretical framework being applicable to any other multiaxial spectral method that needs angular rotations to identify the critical plane.