Practical estimation method for extreme value distribution of von Mises stress in ship structure

This study presents a practical method for estimating the extreme value distribution of von Mises stress for ship structural strength evaluations. The previous methods of calculating this distribution require somewhat complicated numerical calculations, such as multi-dimensional integration. In contrast, the proposed method is based on an asymptotic approximation and can be easily calculated in a similar way to the conventional linear statistical prediction. A closed expression was derived in the case of stress components which have non-zero mean value. The formula is derived under approximations that reflect realistic stress conditions when ships are under severe sea states. Through a structural analysis of a whole ship, it was comprehensively verified that the proposed method has sufficiently high accuracy for structural strength evaluation. Furthermore, a parametric analysis was conducted to clarify its limit of applicability.


General
In the following, (t) = G 1 (t), G 2 (t), G 3 (t) T and (t) = H 1 (t), H 2 (t), H 3 (t) T are dummy variables that represent the vector process: (t) Time derivative of (t) Realization of (t) Mean value of (t) Standard deviation of (t) Covariance matrix between (t) and (t) Correlation coefficients between G i (t) and H j (t) f ( ) Joint probability density function between G 1 (t), G 2 (t), andG 3

(t)
Others c ij Nondimensional coefficient related to the curvature of the isosurface of f ( ) g Gravity acceleration Q Z (z) Extreme value distribution (probability that the maxima of Z(t) exceeds z) T ze Encountered mean zero-upcrossing wave period Upcrossing rate at level Z(t) = z (t) Gaussian vector process of plane stress consisting of normal stress components ( X 1 , X 2 ) and share stress component ( Response amplitude operator of X i (t) (t) Gaussian vector process obtained by linear transformation of (t) y P 1 , y P 2 Coordinates of y 1 and y 2 on plane-P satisfying y 3 = Y 3 y M 1 (z), y M 2 (z) Coordinates of y P 1 and y P 2 in which takes maxima on the isoline Z(t) = z Z(t) Square of von Mises stress Z 0 Value of Z(t) in pre-load z Threshold of Z(t) Wave angle ( = 0 : following sea, = : head sea) Wave spectrum Longitude in -space y P 1 , y P 2 Exponent of f ( ) on plane-P 1 3 Wave frequency e Encounter wave frequency

Introduction
A structural strength evaluation of a ship requires an estimate of the maximum response over the entire period of the ship's operation. The short/long-term prediction method based on linear theory proposed by Fukuda [1] is widely used for ship design because it can be easily calculated from the response amplitude operation (RAO) in the frequency domain and gives reasonable results. However, the method cannot deal with nonlinear quantities such as von Mises stress because it assumes that the response is linear and its extreme value distribution follows a Rayleigh distribution. Even though an estimation of the maximum expected value for von Mises stress is quite important for structural strength evaluations, there are still no methods that are widely adapted for ship structural design. The estimation of the extreme value distribution of von Mises stress can be interpreted as an outcrossing problem on a hypersurface in stress component space, and it is defined by the integral of the probability density function on the hypersurface [2]. The integral cannot be solved analytically in the general case, and numerical approaches are needed to some extent. Although the methods presented in previous studies provide reasonable results [2][3][4], their computational procedures are somewhat complicated, e.g., involving a multi-dimensional numerical integration, and it seems difficult to apply them in practice. For a method to be practical, it must be not only accurate, but also computationally efficient and robust. These features are even more important for ship design because the contributions of all possible sea states must be taken into account when calculating the long-term probability of exceedance.
In terms of computational cost and robustness, the asymptotic approximation of the integral is a very effective approach. There have been several studies on asymptotic formulae of the extreme value distribution of non-Gaussian processes [5][6][7]. However, closed-form expressions for them are limited to special cases. Specifically, a formula for von Mises stress has not been presented for the case that the stress components have non-zero mean value, i.e., when stress in stillwater conditions is considered. Since the still-water loads acting on a ship can be as significant as wave-induced loads, it is essential to consider the mean value of stress components in a structural strength evaluation of a ship.
In light of the above background, the purpose of this paper is to develop a practical method for calculating the extreme value distribution of von Mises stress that is accurate, computationally efficient, and robust. The method assumes that there are three stress components (plane stress conditions) which follow a narrow-band stationary Gaussian process. Here, the author has developed a closed formula for the extreme value distribution that is based on an asymptotic expansion of the integral by an extension of Laplace's method. In particular, the proposed formula takes into account the non-zero mean values in closed form, which has not been shown in previous studies [5][6][7]. In deriving the formula, approximations are applied to the extent that accuracy does not deteriorate under realistic stress conditions in which the ship is under severe sea states. The applicability of proposed method is verified comprehensively by conducting a whole ship structural analysis in waves.

Formulation of extreme value distribution of von Mises stress
First, we define the random process of stress components and von Mises stress. Next, for the sake of simplicity, we introduce a variable transformation of the stress components and express von Mises stress in a sum-of-squares form. Then, we define the extreme value distribution in terms of the upcrossing rate and finally give a specific formulation of the extreme value distribution of von Mises stress using stochastic parameters of the stress components.

Definition of a random process
Let X 1 (t) and X 2 (t) be normal stress components and X 3 (t) be the shear stress component of the plane stress at a certain position in a ship structure. Furthermore, let us assume them to be random variables following a stationary Gaussian distribution in short-term sea states. Their mean values X i ∶= E X i (t) are interpreted as the stress occurring in the still-water condition. Voigt's notation for them, (t) ∶= X 1 (t), X 2 (t), X 3 (t) T , is a Gaussian vector process with mean value vector ∶= X 1 , X 2 , X 3 T and covari- The joint probability density function (PDF) of (t) is The state is also denoted as (t) ∼ N , . The time derivative of (t) , denoted as ̇ (t) , is a stationary Gaussian vector process with zero mean, i.e., ̇ (t) ∼N 0, ̇ ̇ . In addition, the covariance matrix between (t) and ̇ (t) is Since the Gaussian process X i (t) and its time derivative Ẋ i (t) are mutually independent, the diagonal components of ̇ are all zero.
The square of the von Mises stress Z(t) is expressed in the following quadratic form of the stress components.
The symbol for the time variable (t) of the random process is omitted hereafter.

Variable transformation
To simplify the extreme value distribution of Z in Eq. (2), we will transform the stress component vector into another Gaussian vector process such that Z can be expressed in sum-of-squares form, according to the procedure presented by Segalman [8,9].
The vector process follows ∼ N , , where is the following diagonal matrix: The variable transformation → can be applied when Z is a quadratic form of . Appendix 1 shows the procedure of the transformation, which is a modification of Segalman's in which only one diagonalization is required. Furthermore, Appendix 1 shows how to obtain the covariance matrixes , ̇ , and ̇ from the RAOs of the stress components for the case of a ship has advancing speed in short-crested irregular waves.
The mean value of Z can be expressed as follows.
where Z 0 means the value of Z in pre-load, i.e., in still-water condition. (2)

Definition of extreme value distribution
As it is well known, the number of maxima can be approximated by the number of upcrossings for a narrow-band process whose minima are not so high in value. Therefore, the extreme value distribution Q Z (z) , i.e., the probability that an extreme value exceeds Z = z , can be approximated as follows.
where T p is the mean period between maxima of Z , and v + Z (z) is the upcrossing rate at level z . The definition of (2.8), however, is not suitable for a statistical prediction of von Mises stress, because it is difficult to formulate T p . In addition, the number of maxima of von Mises stress is about twice that of the maxima of the stress component; hence, definition by Eq. (8) is a different index from the exceedance probability of the conventional short/long-term prediction method.
Instead, the following definition is used in this study.
where T ze is the encountered mean zero-upcrossing wave period, which considers the effect of the ship advancing into waves, where is the encounter wave frequency and Φ , , , U , and g are the wave spectrum, wave frequency, wave angle, ship speed, and gravity acceleration, respectively. T ze is approximately equal to the mean zero-upcrossing period of the fluctuation components of stress X i − X i . As defined in Eq. (9), the maximum expected value of Z occurring once in N-waves coincides with z satisfying Q Z (z) = 1∕N and can be treated in the same way as the index of the exceedance probability used in the conventional short/ long-term prediction method. Strictly speaking, it is not accurate to call it an exceedance probability because Eq. (9) is an expected value of upcrossings in terms of T ze and can exceed 1. However, for the purpose of estimating the maximum expected value of Z in a certain duration, which is used for short/long-term prediction, definition by Eq. (9) is appropriate and the fact that Q z (z) exceeds 1 does not cause any problems.

Formulation of the upcrossing rate
Upcrossing rate v + Z (z) can be calculated using Rice's formula [10]: where f ZŻ (z,̇z) is the joint PDF of Z and its time derivative Ż . Gupta proposed a direct way to find f ZŻ (z,̇z) for von Mises stress [4].
Veneziano generalized Rice's formula (12) and presented the following equation for the case, where Z is expressed as a general scalar function of X i [2].
in the direction normal to the ellipsoid and directed outward, and ḟX n | = ̇x n | is the conditional PDF of Ẋ n . Madsen derived a specific form for Eq. (13) in the case that Z is the square of von Mises stress and consists of three stress components [3]. For Madsen's expression, using the variable in Eq. (4) instead of , the upcrossing rate can be expressed as Although it is an exact formula without any special assumptions, this formula is not so practical because it requires a numerical solution of the double integral. Moreover, it is difficult to ensure its reliability for cases where the matrix ) is close to being singular, as is often the case.

Asymptotic formula of the upcrossing rate
This section derives a practical closed formula for the upcrossing rate v + Z (z) by adapting an asymptotic expansion and making some approximations of the integral in Eq. (14).

Formulation assuming that components are mutually independent
Let us assume that Y and their time derivatives ̇ are mutually independent. In this case, and Eq. (14) can be simplified as where It can be confirmed that the same result is derived by Fukuda's derivation procedure [11], which makes the same approximations as in Eqs. (20) and (21). It has been theoretically shown by Hagen that the effect of the approximation in Eq. (20) on v + Z (z) is secondary, and the loss of accuracy can be ignored [6]. On the other hand, Eq. (21) is not general and does not hold in all cases. However, as we will see in Sect. 4, Eq. (21) holds in most cases without any problem, because when is orthogonally transformed, the covariance matrix ̇ ̇ is almost diagonalized.
When Y 2 ∕ Y 1 and Y 3 ∕ Y 1 are close to zero, the exponential part of Eq. (22) changes sharply, making it difficult to ensure the accuracy of the numerical integration. In particular, Y 3 ∕ Y 1 often takes a very small value. Consequently, we next derive the asymptotic formulae of Eq. (22) for the case where Y 3 ∕ Y 1 are close to zero.
When Y 3 ∕ Y 1 → 0 , the integral on the right-hand side of Eq. (22) is dominated by the contribution of the integrand from the vicinity of the plane-P formed by y 3 = Y 3 . In this case, except for the exponential function of y 3 , the change with respect to y 3 is gradual, so it can be represented by the value on the plane-P. Hence, the following asymptotic approximation is derived by applying the Gaussian integral formula.
Accordingly, the upcrossing rate asymptotically follows as where y P 1 and y P 2 are the coordinates of y 1 and y 2 on plane-P, as shown in Fig. 1.
. Expression of Eq. (30) is simple and is suitable for a rough order evaluation of z . When we set Q Z (z) ≅ T zY 1 v + Z (z) and solve for z , we get This formula is not accurate, but it is useful as a reference value for calculating the distribution of Q Z (z) while varying z.

Asymptotic expansion of Eq. (25)
Here, we derive a more practical formula without an integral or case separation on the basis of an asymptotic expansion using Laplace's method.
Even though the Laplace's method is a series representation of the Laplace integral under the assumption that z is sufficiently large, it is known to be a good approximation over a wide range even if only the first term is considered [14]. When = 0 , the integral in Eq. (22) becomes the Laplace integral, so we can directly apply Laplace's method. Appendix 2 derives the following asymptotic formula: is a nondimensional coefficient related to the curvature of the isosurface of f ( ) on y i y j -plane. Equation (33) coincides with Breitung's result [5].
On the other hand, the derivation is not straightforward when ≠ 0 , because it is difficult to explicitly express the coordinates of the reference points for the Taylor expansion. In this study, the coordinates of the reference points are approximated. Equation (25) is used as the base formula of the asymptotic expansion, because Y 3 ∕ Y 1 usually does not take a large value, as will be shown in Sect. 4.
In taking the asymptotic expansion of Eq. (25), we must find the angle M where is maximized on the isosurface of Z (the circle formed by the intersection of the plane-P and the isosurface of Z ) and consider the contribution of the integrand only in the vicinity of M . In the case of the integral in Eq. (25), there are usually two maxima of and the contribution of the second largest maximum cannot be ignored in the case Y 1 ≪ Y 1 . Therefore, to consider both two maxima, let us denote the (y P 1 , y P 2 )-coordinates where is at a maximum by y M+ 1 , y M+ where The coefficient √ c 31 is additionally taken into account to consider the effect of Y 3 ∕ Y 1 by referring to the asymptotic formula (33)with zero mean.
Next, we find the points (y M± 1 , y M± 2 ) where reaches maximal value. They are at the intersections of the isoline of z and the curve on which gradient of (y P 1 , y P 2 ) is directed to the origin of the plane-P, that is, Equation (40) is the hyperbola that passes through the origin y P 1 , y P 2 = (0, 0) as well as the mean value point and is asymptotic to the two lines: y P 1 = c 12 Y 1 , and y P 2 = c 21 Y 2 . The schema of the contour of z and on plane-P and hyperbola of Eq. (40) which takes an extreme value is shown in Fig. 2. From Fig. 2, it can be understood that has two maxima and two minima on the contour line of Z and the maximum is in the region of sign( Y 1 ) y P 1 ≥ 0. To obtain an explicit expression for (y M± 1 , y M± 2 ) , which is defined as the points that satisfy Eqs. (38) and (40) simultaneously, it is necessary to solve a quartic function and the closed formula of (y M± 1 , y M± 2 ) becomes quite complicated. Therefore, we will treat this problem approximately. First, for simplicity, let us approximate the hyperbola of Eq. (40) as a symmetric curve which mirrors the hyperbola on the side of sign( Y 1 ) y P 1 ≥ 0 with respect to the y P 2 -axis, i.e., y M+ Then, as an explicit function of z , y M 2 can be approximated as follows: This formula has the following value when Y 1 = 0 , as can be seen by taking the limit Y 1 → 0.
The derivation of Eq. (41) is shown in Appendix 3. The approximation formula (41) is a curve whose asymptotic behavior as z → ∞ is correct on the order of O(1∕ ) compared to the true solution, and the curve passes through the origin and the mean value point ( Y 1 , Y 2 ) at z = 0 and z = Z 0 , respectively. Furthermore, it always satisfies y M 2 < when z ≥ Z 0 so that the square root of y extreme value distribution Q Z (z) , Eq. (46) requires only Y i and Y i when Q Z (z) ≅ T zY 1 v + Z (z) holds (in Sect. 4 indicates that this approximation is usually valid).

Target ship structural model for analysis
To validate the applicability of the proposed extreme value distribution formula (46) to ship structural design, a direct load and structure analysis of a whole ship FEM model was conducted using the DLSA-Basic system [15], wherein wave-induced stresses and still-water stresses were obtained of all shell elements. The target model was a bulk carrier with L = 280m, and the full load condition with homogeneous loading (homo.) and alternative loading (alt.) were considered. For each condition, conventional short/long-term predictions were performed for the three stress components in the neutral plane of all shell elements, and the top 100 elements of the long-term maximum expected value (including the still-water component) were selected for each sevenelement group where the mean wave angle in the MSS (Most Severe Short-term sea states) [16] was 0°, 30°, …, 180°. The selected 700 elements can be regarded as the primary structural members affected by various load factors in all wave directions and are considered suitable for comprehensive verification. Figure 3 shows the ship FEM model with the selected elements highlighted in the homo./alt. conditions.

Histogram of stochastic parameters of selected elements in the most severe short-term sea state
Before the numerical validation, in order to capture a realistic range of stochastic parameters of the stress vector process, the mean value and covariance matrixes , ̇ , ̇ ̇ of the vector process in the MSS, which is converted from the stress component vector as presented in Appendix 1, were calculated for each of the 2 loading conditions × 700 elements. Figure 4 shows scatter plots on the Y 2 ∕ Y 1 -Y 3 ∕ Y 1 plane of 700 elements × 2 conditions. From the figure, almost all values of Y 2 ∕ Y 1 are less than 0.5 and almost all values of Y 3 ∕ Y 1 are less than 0.1 (maximum 0.33). When the standard deviation ratio Y 2 ∕ Y 1 is close to 1, the element is affected by two or more uncorrelated and comparable dominant load factors, e.g., vertical bending moment and roll motion under the MSS. The fact that Y 3 ∕ Y 1 is generally very small implies that the hull structural members are dominated by only two uncorrelated (and comparable) load factors. This fact might also hold in the case of a 3-dimensional stress condition. That is, when the six standard deviations Y 1 , … , Y 6 are obtained from the orthogonal transformation of six stress components of solid element and arranged in descending order, Y i ∕ Y 1 ( i ≥ 3 ) might be negligible. Hence, we infer that the proposed formula (46), in which is defined ) , is also applicable to a solid element. Figure 5 shows histograms of ratios of the mean value of Y i against the square root of the mean value of Z . Since the mean value of Z is Z =

as written in
Eq. (6), the fact that the ratio is close to ±1 means that the mean value Y i is very large compared to the wave-induced component. Figure 5 confirms that their ratios are in wide range and are close to ±1 for some elements; therefore, it is essential to consider the mean value to estimate the maximum value of stress in the ship structure. These absolute values are roughly in the order , which is due to the correlation between the wave-induced load and still-water load. In addition, the absolute value of the alt.condition is larger than that of the homo.-condition because under the alternative loading, high stress is likely to occur even in the still-water condition. Figure 6 shows histograms of the ratios of the mean zeroupcrossing period between Y i -component T zY i ( ∶= 2 Y i ∕̇Y i ) and incident wave T ez . Here, it can be seen that T zY i takes a value close to T ez and differs by 20% at most. Therefore, Fig. 4 Scatter plots on plane of the standard deviation ratio Y 2 ∕ Y 1 and Y 3 ∕ Y 1 for the 700 selected elements (left: homo., right: alt.)

Fig. 5 Histograms of mean values Y i normalized by square root of
) ) for the 700 selected elements (left: homo., right: alt.) considering that T ez in Q Z (z) has a logarithmic effect on z , it can be approximated as T ez ≅ T zY i . Furthermore, we can say that The reason why the period of the main component T zY 1 tends to be longer than T ez is that the dominant load factor has a slightly longer period than the wave period in the MSS. This depends on the scale of the ship, and T zY 1 is expected to be shorter than T ez for smaller ships.

Figures 7 and 8 show histograms of correlation coefficients between
. It is found from Fig. 7 that the absolute value of Ẏ iẎj is at most 0.3; hence, Eq. (21) is a reasonable approximation. On the other hand, Fig. 8 shows that Eq. (20) does not hold because Y iẎj can be large. However, the effect of Y iẎj is secondary for the upcrossing rate [6] and thus Eq. (20) does not cause any problems, as will be confirmed numerically below.

Generation of time series for von Mises Stress
The extreme value distribution obtained from the generated time series of von Mises stress was also used for verification. The time series of were generated from the RAOs of the stress components X i and wave spectrum Φ ( , ) with directional scatter, as follows. (47) where en ∶= n − 2 n g U cos m , and n,m , the phase advance of the incident wave component, is a uniform random number in the range [0, 2 ] . The time series of the von Mises stress was calculated from the time series of the stress components at each time step, and its extreme value distribution was obtained by picking up the maxima. The ISSC wave spectrum was used for Φ ( , ) [17]. Assuming an event in a short-term sea state, the duration of the time series and the number of repetitions are sufficiently large to numerically stabilize the extreme value distribution. Specifically, the spectrum was divided into 1000 equal areas against wave frequency ( N = 1000 ) and into 5 parts against wave direction ( M = 5 ). The duration of the time series was 2000 T ez [s], and it was repeated 50 times with different random numbers of n,m .

Comparison of extreme value distributions of von mises stress
We compared the proposed formula (46) for the extreme value distribution of Mises stress with the Madsen's exact formula (14) and the extreme value distribution obtained The plot of the extreme value distribution Q Z (z) for each method is shown in Figs. 9, Fig. 10, and 11 for elements A, B, and C, respectively. The range of the horizontal axis is limited to Z ≥ Z 0 , because Z does not have to be estimated in the region smaller than Z 0 for practical use, and the upcrossing rate v + Z (z) peaks near Z = Z 0 and is different from the actual distribution of Q Z (z) in Z < Z 0 . From these figures, it can be seen that the proposed formula (46) is in good agreement with the Madsen's formula and the results of the time series. The wobble in the time series distribution for Q < 1∕1000 is due to statistical instability. The shape of the distribution of element B differs from those of elements A and C. This can be understood from the following rough estimation formula obtained by Eq. (32): Figs. 12 and 13 compare Madsen's formula and the proposed formula for the level z at which Q Z (z) = 1∕10 and 1∕1000 for the two conditions and 700 selected elements. Since the distribution of extremes in the time series showed variations due to randomness, the results are compared with those of Madsen's formula (14), which is completely specified by the response spectrum. From these figures, it can be seen that the proposed formula (46) has very high accuracy. The small error at Q Z (z) = 1∕10 was mainly caused by the asymptotic expansion, while the approximation by Eqs. (20) and (21) had extremely small error. In other words, there is negligible difference between the Madsen's formula (14) and the integral on the plane-P in Eq. (22).

Applicable range of proposed formula
Finally, we investigated the applicable parameter range of the proposed formula (46). Here, we investigated the error distribution of the proposed formula in a five-parameter space: standard deviation ratios The error rate is defined as where z Formula Q=1∕1000 is the value of z obtained from the proposed formula (46) that satisfies Q Z (z) = 1∕1000 , and z Exact Q=1∕1000 is the value obtained from Madsen's formula (14). Figure 14 shows the contours of the error rate on the  Fig. 11 Comparison of exceedance probabilities of peak value of Z at element C. in Fig. 14 is set to 0.99 because the proposed formula goes to infinity as Y 2 ∕ Y 1 → 1 , but this singularity does not cause a practical problem. Figure 15 shows the error rate on the Figure 14 and 15 indicate that the accuracy of 2% is guaranteed for a wide range of parameters. However, in Fig. 14, when Y 1 = 0 and Y 2 = 3 Y 3 (upper right and lower right figures), the proposed formula over-estimates in the region distance between the two points ( y M± 1 , y M± 2 ), where takes maximal value as shown in Fig. 2, becomes shorter, and these contributions overlap in the asymptotic Eq. (35). However, it is found from Fig. 15 that the region where has a large value is very close to the Y 1 = 0 axis and it is extremely rare for both Y 1 ≅ 0 and Y 2 ∕ Y 1 > 0.85 to be true (see Figs. 4 and 5). Therefore, although the proposed formula over-estimates by more than 5% under the condition  Fig. 15) and Y 2 ∕ Y 1 > 0.85 , there is little problem in practice.

Conclusion
The author proposed a practical method for estimating the extreme value distribution of von Mises stress that can be easily applied to ship design. The proposed method combines the variable transformation of stress components presented by Segalman and an asymptotic approximation for the integral of the upcrossing rate. The features of the proposed method are enumerated as follows.
(i) A closed formula for the upcrossing rate is given by making an asymptotic approximation of the integral in the case of stress components with non-zero mean value. The formula is based on an idea like that of Laplace's method, which is extended to the non-zero mean value case and the solution of the quartic equation of the reference point is simplified for an asymptotic expansion. Other approximations within realistic parameter ranges are applied; e.g., the transformed vector process and its time derivative are mutually independent and the standard deviation ratio between the variables is not large. (ii) The computational procedure of the proposed method does not differ much from the conventional linear short/ long-term statistical prediction. The only additional calculations are for deriving the covariance matrix of the stress components from the RAOs and wave spectrum and diagonalizing the 3 × 3 covariance matrix. The proposed method does not require the numerical integration commonly required by other methods; thus, robust results can be obtained at a low computational cost. (iii) Through the structural analysis of a whole structural model of a bulk carrier in waves, it was confirmed that the proposed method has sufficient accuracy for the maximum expected value of von Mises stress. Furthermore, the parametric study clarified the application range of the proposed method and confirmed that the accuracy would not deteriorate for structures in waves.

Appendix 1. Variable Transformation of von Mises stress into Some of Squares Form
This appendix describes the procedure to obtain the covariance matrixes , ̇ ̇ , and ̇ from the RAOs of and the wave spectrum of short-crested irregular waves. Using these matrices, is transformed into another vector process in which Z is expressed as a sum of squares.
Calculation of covariance matrix of X Let us denote the RAO of X i as X i ( , ) (complex number), where and are wave frequency and wave angle, respectively, and denote the wave spectrum of short-crested irregular waves as Φ ( , ) . The components of the covariance matrix in the shortterm sea state can be calculated as where the superscript "*" denotes the complex conjugate, Similarly, ̇ ̇ and ̇ can be derived by replacing X → i eX in Eq. (51). The above equation defines the wave angle such that = 0 is a following sea and = is a head sea. Thus, ̇ ̇ and ̇ are calculated as

Sum of squares expression of von Mises stress
The procedure of standard normalization of the joint PDF of vector process is commonly applied, but in addition, it is also possible to make its derivatives ̇ independent by orthogonal transformation. That is, without loss of generality, vector process can be transformed into another vector process which satisfies = and ̇ ̇ = diag( 2V ) as can be seen in, e.g., Ref. [2,6,7]. However, in this study, we follow Segalman's method [8,9], which standardizes the iso-ellipsoid of Z instead of the PDF of ̇ . The procedure presented in this study is a modification of Segalman's method in which only one diagonalization is required. First, we transform so that the iso-ellipsoid of Z in -space becomes a sphere, as follows.
where (51) Σ X i X j = ∫ − ∫ is an orthogonal matrix which diagonalizes , and is a diagonal matrix where A i are the eigenvalues of , i.e., = T . Since Z in Eq. (2) can be expressed by a vector process , the isosurface of Z is a sphere in -space. The covariance matrix of can be obtained from as can be understood from the transformation of the exponent portion of the joint PDF of , Eq. (1).
Next, we transform into a vector process whose components are mutually independent by performing an orthogonal transformation. To do so, the diagonalization was conducted using an orthogonal matrix as follows.
Here, the order is Y 1 ≥ Y 2 ≥ Y 3 . can be obtained by transforming as follows.
The vector process follows ∼ N , , where = T T and is a diagonal matrix as shown in Eq. (58). Furthermore, Z is expressed as a sum of squares of Y i : The covariance matrix of ̇ , that is ̇ ̇ , can also be derived by transforming ̇ ̇ as The same transform is applicable to ̇ and ̇ . Furthermore, the mean value of Z is expressed as the standard deviation and mean value of , i.e., Thus, it turns out that a nonlinear quantity Z expressed in any quadratic form can be expressed as a sum of squares of mutually independent Gaussian vector process which .
(59) = T = T T . (60) (61) ̇ ̇ = T T ̇ ̇ . (62) is obtained by making a linear transformation of . A schematic diagram of the above derivation process for the two-dimensional case is shown in Fig. 16. Here, it can be understood that this method is applicable when the isosurface of Z is an ellipsoid, i.e., that Z is in quadratic form.

Appendix 2. Asymptotic expansion by Laplace's method
General case Let us consider the bi-variable Laplace integral, where f ( , ) and g(z, , ) are continuous real functions and g(z, , ) is an increasing function with respect to z . When g(z, , ) takes a maximum value at a point ( , ) = 0 , 0 in the integration range, it can be asymptotically expanded using Laplace's method [14]. In this case, when z is sufficiently large, the integrand is dominated by the contribution near ( , ) = 0 , 0 . Therefore, substituting the Taylor expansion of g(z, , ) around ( , ) = 0 , 0 , Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.