# An Analytical Solution for Predicting the Vibration-Fatigue-Life in Bimodal Random Processes

Authors: Chaoshuai Han; Yongliang Ma; Xianqiang Qu; Mindong Yang
Journal: Shock and Vibration
 (2017)

---

## Abstract

Predicting the vibration-fatigue-life of engineering structures subjected to random loading is a critical issue for. Frequency methods are generally adopted to deal with this problem. This paper focuses on bimodal spectra methods, including Jiao-Moan method, Fu-Cebon method, and Modified Fu-Cebon method. It has been proven that these three methods can give acceptable fatigue damage results. However, these three bimodal methods do not have analytical solutions. Jiao-Moan method uses an approximate solution, Fu-Cebon method, and Modified Fu-Cebon method needed to be calculated by numerical integration which is obviously not convenient in engineering application. Thus, an analytical solution for predicting the vibration-fatigue-life in bimodal spectra is developed. The accuracy of the analytical solution is compared with numerical integration. The results show that a very good agreement between an analytical solution and numerical integration can be obtained. Finally, case study in offshore structures is conducted and a bandwidth correction factor is computed through using the proposed analytical solution.


---

## Body

## 1. Introduction

Engineering structures from different fields (e.g., aircrafts, wind energy utilizations, and automobiles) are commonly subjected to random vibration loading. These loads often cause structural fatigue failure. Thus, it is significant to carry out a study on assessing the vibration-fatigue-life [1, 2].Vibration fatigue analysis commonly consists of two part processes: structural dynamic analysis and results postprocessing. Structural dynamic analysis provides an accurate prediction of the stress responses of fatigue hot-spots. Once the stress responses are obtained, vibration fatigue can be successfully performed. Existing technologies such as operational modal analysis [3], finite element modeling (FEM), and accelerated-vibration-tests are mature and applicable to obtain the stress of structures [4, 5]. Therefore, the crucial part of a vibration fatigue analysis focuses on results postprocessing.The postprocessing is usually used to calculate fatigue damage based on known stress responses. When the stress responses are time series, fatigue can be evaluated using a traditional time domain method. However, the stress responses of real structures are mostly characterized by the power spectral density (PSD) function. Thus, frequency domain method becomes popular in vibration fatigue analysis [6, 7].A bimodal spectrum is a particular PSD in the random vibration stress response of a structure. For some simple structures, the stress response of structures will show explicit characterization of two peaks. One peak of the bimodal spectrum is governed by the first-order natural frequency of the structure; another is dominated by the main frequency of applied loads. Therefore, some bimodal methods for fatigue analysis can be adopted [8–10]. Moreover, several experiments (e.g., vibration tests on mechanical components) and numerical studies (e.g., virtual simulation of dynamic using FEM) also obtain have shown that the stress PSD is a typical bimodal spectrum [5, 7, 11]. However, for some complex flexible structures, the PSD of the stress response of structures usually is a multimodal and wide-band spectrum. For this situation, existing general wide-band spectral methods such as Dirlik method [12], Benasciutti-Tovo method [10], and Park method [13, 14] can be used to evaluate the vibration-fatigue-life. Recently, Braccesi et al. [15, 16] proposed a bands method to estimate the fatigue damage of a wide-band random process in the frequency domain. In order to speed up the frequency domain method, Braccesi et al. [17, 18] developed a modal approach for fatigue damage evaluation of flexible components by FEM.For fatigue evaluation in bimodal processes, some specific formulae have been proposed. Jiao and Moan [8] provided a bandwidth correction factor from a probabilistic point of view, and the factor is an approximate solution derived by the original model. The approximation inevitably leads to some errors in certain cases. Based on an similar idea, Fu and Cebon [9] developed a formula for predicting the fatigue life in bimodal random processes. In the formula, there is a convolution integral. The author claimed that there is no analytical solution for the convolution integral which has to be calculated by numerical integration. Benasciutti and Tovo [10] compared the above two methods and established a Modified Fu-Cebon method. The new formula improves the damage estimation, but it still needs to calculate numerical integration.In engineering application, the designers prefer an analytical solution rather than numerical methods. Therefore, the purpose of this paper is to develop an analytical solution to predict the vibration-fatigue-life in bimodal spectra.


## 2. Theory of Fatigue Analysis

### 2.1. Fatigue Analysis

The basicS-N curve for fatigue analysis can be given as(1)N=K·S-m,where S represents the stress amplitude; N is the number of cycles to fatigue failure, and K and m are the fatigue strength coefficient and fatigue strength exponent, respectively.The total fatigue damage can then be calculated as a linear accumulation rule after Miner [19](2)D=∑niNi=1K∑niSim,where ni is the number of cycles in the stress amplitude Si, resulting from rainflow counting (RFC) [20], and Ni is the number of cycles corresponding to fatigue failure at the same stress amplitude.When the stress amplitude is a continuum function and its probability density function (PDF) isfS(S), the fatigue damage in the duration T can be denoted as follows:(3)D=νc·TK∫0∞Sm·fSSdS,where νc is the frequency of rainflow cycles.For an ideal narrowband process,fS(S) can be approximated by the Rayleigh distribution [21]; the analytical expression is given as (4)fSS=Sλ0exp⁡-S22λ0.Furthermore, the frequency of rainflow cycles νc can be replaced by rate of mean zero upcrossing ν0.According to (3), an analytical solution of fatigue damage [22] for an ideal narrowband process can be written as(5)DNB=ν0·TK2λ0mΓm2+1,where Γ· is the Gamma function.For general wide-band stress process, fatigue damage can be calculated by a narrowband approximation (i.e., (5)) first, and bandwidth correction is made based on the following model [23]:(6)DWB=ρ·DNB.In general, bimodal process is a wide-band process; thus, the fatigue damage in bimodal process can be calculated through (6).


### 2.2. Basic Principle of Bimodal Spectrum Process

Assume that a bimodal stress processX(t) is composed of a low frequency process (LF) XL(t) and a high frequency process (HF) XH(t). (7)Xt=XLt+XHt,where XL(t) and XH(t) are independent and narrow Gaussian process.The one-sided spectral density function ofX(t) can be summed from the PSD of LF and HF process. (8)Sω=SLω+SHω.Theith-order spectral moments of S(ω) are defined as(9)λi=∫0∞ωi·SLω+SHωdω=λi,L+λi,H.The rate of mean zero upcrossing corresponding toXL(t) and XH(t) is(10)ν0,L=12πλ2,Lλ0,L,ν0,H=12πλ2,Hλ0,H.The rate of mean zero upcrossing ofX(t) can be expressed as(11)ν0=12πλ2,L+λ2,Hλ0,L+λ0,H=ν0,L2λ0,L+ν0,H2λ0,Hλ0,L+λ0,H.According to (5), (9), and (11), narrowband approximation of bimodal stress process X(t) can be given as(12)DNB,X=ν0,L2λ0,L+ν0,H2λ0,Hλ0,L+λ0,HTK2λ0,L+2λ0,HmΓm2+1.Equation (12) is known as the combined spectrum method in API specifications [24].The existing bimodal methods proposed by Jiao and Moan, Fu and Cebon, and Benasciutti and Tovo are based on the idea: two types of cycles can be extracted from the rainflow counting, one is the large stress cycle, and the other is the small cycle [8–10]. The fatigue damage due to X(t) can be approximated with the sum of two individual contributions.(13)D=Dl+Ds,where Dl represents the damage due to the large stress cycle and Ds denotes the damage due to the small stress cycle.


## 2.1. Fatigue Analysis

The basicS-N curve for fatigue analysis can be given as(1)N=K·S-m,where S represents the stress amplitude; N is the number of cycles to fatigue failure, and K and m are the fatigue strength coefficient and fatigue strength exponent, respectively.The total fatigue damage can then be calculated as a linear accumulation rule after Miner [19](2)D=∑niNi=1K∑niSim,where ni is the number of cycles in the stress amplitude Si, resulting from rainflow counting (RFC) [20], and Ni is the number of cycles corresponding to fatigue failure at the same stress amplitude.When the stress amplitude is a continuum function and its probability density function (PDF) isfS(S), the fatigue damage in the duration T can be denoted as follows:(3)D=νc·TK∫0∞Sm·fSSdS,where νc is the frequency of rainflow cycles.For an ideal narrowband process,fS(S) can be approximated by the Rayleigh distribution [21]; the analytical expression is given as (4)fSS=Sλ0exp⁡-S22λ0.Furthermore, the frequency of rainflow cycles νc can be replaced by rate of mean zero upcrossing ν0.According to (3), an analytical solution of fatigue damage [22] for an ideal narrowband process can be written as(5)DNB=ν0·TK2λ0mΓm2+1,where Γ· is the Gamma function.For general wide-band stress process, fatigue damage can be calculated by a narrowband approximation (i.e., (5)) first, and bandwidth correction is made based on the following model [23]:(6)DWB=ρ·DNB.In general, bimodal process is a wide-band process; thus, the fatigue damage in bimodal process can be calculated through (6).


## 2.2. Basic Principle of Bimodal Spectrum Process

Assume that a bimodal stress processX(t) is composed of a low frequency process (LF) XL(t) and a high frequency process (HF) XH(t). (7)Xt=XLt+XHt,where XL(t) and XH(t) are independent and narrow Gaussian process.The one-sided spectral density function ofX(t) can be summed from the PSD of LF and HF process. (8)Sω=SLω+SHω.Theith-order spectral moments of S(ω) are defined as(9)λi=∫0∞ωi·SLω+SHωdω=λi,L+λi,H.The rate of mean zero upcrossing corresponding toXL(t) and XH(t) is(10)ν0,L=12πλ2,Lλ0,L,ν0,H=12πλ2,Hλ0,H.The rate of mean zero upcrossing ofX(t) can be expressed as(11)ν0=12πλ2,L+λ2,Hλ0,L+λ0,H=ν0,L2λ0,L+ν0,H2λ0,Hλ0,L+λ0,H.According to (5), (9), and (11), narrowband approximation of bimodal stress process X(t) can be given as(12)DNB,X=ν0,L2λ0,L+ν0,H2λ0,Hλ0,L+λ0,HTK2λ0,L+2λ0,HmΓm2+1.Equation (12) is known as the combined spectrum method in API specifications [24].The existing bimodal methods proposed by Jiao and Moan, Fu and Cebon, and Benasciutti and Tovo are based on the idea: two types of cycles can be extracted from the rainflow counting, one is the large stress cycle, and the other is the small cycle [8–10]. The fatigue damage due to X(t) can be approximated with the sum of two individual contributions.(13)D=Dl+Ds,where Dl represents the damage due to the large stress cycle and Ds denotes the damage due to the small stress cycle.


## 3. A Review of Bimodal Methods

### 3.1. Jiao-Moan (JM) Method

To simplify the study,X(t), XL(t), and XH(t) are normalized as X∗(t), XL∗(t), and XH∗(t) through the following transformation:(14)X∗t=Xtλ0=XLtλ0+XHtλ0=XL∗t+XH∗tand then(15)λ0∗=λ0,L∗+λ0,H∗=1,where(16)λ0,L∗=λ0,Lλ0,λ0,H∗=λ0,Hλ0.Jiao-Moan points out that the small stress cycles are produced by the envelope of the HF process, which follows the Rayleigh distribution. The fatigue damage due to the small stress cycles can be obtained according to (5).While the large stress cycles are from the envelop process,P(t) (see Figure 1), the amplitude of P(t) is equal to(17)Qt=RLt+RHt,where RL(t) and RH(t) are the envelopes of XL∗(t) and XH∗(t), respectively.Figure 1
Bimodal processX∗(t), the envelope process P(t), and the amplitude process Q(t).The distribution ofQ(t) can be written as a form of a convolution integral(18)fQq=∫0qfRLq-xfRHxdx=∫0qfRLyfRHq-ydy.RL(t) and RH(t) obey the Rayleigh distribution; therefore, (18) has an analytical solution which is given [8](19)fQq=qλ0,L∗·exp⁡-q22λ0,L∗+qλ0,H∗·exp⁡-q22λ0,H∗+exp⁡-q22·2πλ0,L∗λ0,H∗·Φqλ0,L∗λ0,H∗+Φqλ0,H∗λ0,L∗-1·q2-1.The rate of mean zero upcrossing due toP(t) can be calculated as(20)ν0,P=λ0,L∗ν0,L1+λ0,H∗λ0,L∗ν0,Hν0,LδH2,where(21)δH=1-λ1,H∗2λ0,H∗λ2,H∗.An approximation was made by Jiao and Moan for (19) as follows [8]:(22)fQq≈λ0,L∗-λ0,L∗λ0,H∗·q·exp⁡-q22λ0,L∗+2πλ0,L∗λ0,H∗·q2-1·exp⁡-q22.After the approximation, a closed-form solution of the bandwidth correction factor can be then derived [8](23)ρ=ν0,Pν0×λ0,L∗m/2+21-λ0,H∗λ0,L∗+πλ0,L∗λ0,H∗mΓm/2+1/2Γm/2+1+ν0,Hν0λ0,H∗m/2.Finally, the fatigue damage can be obtained as (6) and (12).


### 3.2. Fu-Cebon (FC) Method

Similarly to JM method, Fu and Cebon also considered that the total damage is produced by a large cycle (SH+SL) and a small cycle (SL), as depicted in Figure 2. The small cycles are from the HF process, and the distribution of the amplitude PSs(S) is a Rayleigh distribution, as shown in (4). However, the number of cycles associated with the small cycles ns is different from JM method and equals ν0,H-ν0,L·T. According to (5), the damage due to the small cycles is(24)Ds=ν0,H-ν0,L·TK2λ0,HmΓ1+m2.Figure 2
The large cycles and small cycles for a random stress process.The amplitude of the large cyclesSl can be approximated as the sum of amplitude of the LF and HF processes, the distribution of which can be expressed by a convolution of two Rayleigh distributions [9]. (25)PSlS=∫0SPSLyPSHS-ydy=∫0SPSLS-yPSHydy=1λ0,Lλ0,He-S2/2λ0,H∫0SSy-y2e-Uy2+VSydy,where U=1/2λ0,L+1/2λ0,H and V=1/λ0,H.The number of cycles of the large cycles isnl=ν0,L·T. Thus, the fatigue damage due to the large stress cycles can be expressed by(26)Dl=ν0,L·TK∫0∞SmPSlSdS.Equation (26) can be calculated with numerical integration [9, 10]. Therefore, the total damage can be obtained according to (13).


### 3.3. Modify Fu-Cebon (MFC) Method

Benasciutti and Tovo made a comparison between JM method and FC method and concluded that using the envelop process is more suitable [10]. Thus, a hybrid technique is adopted to modify the FC method. More specifically, the large cycles and small cycles are produced according to the idea of FC method. The number of cycles associated with the large cycles is defined similarly to JM method. That is, nl=ν0,P·T, while the number of cycles corresponding to the small cycles is ns=ν0,H-ν0,P·T. The total damage for MFC method can be then written according to (13).Although the accuracy of the MFC method is improved, the fatigue damage still has to be calculated with numerical integral.


### 3.4. Comparison of Three Bimodal Methods

Detailed comparison of the aforementioned three bimodal methods can be found in Table1. In all methods, the amplitude of the small cycle obeys Rayleigh distribution, and the corresponding fatigue damage has an analytical expression as in (5); the distribution of amplitude of the large cycle is convolution integration of two Rayleigh distributions, and the relevant fatigue damage can be calculated by (26).Table 1
Comparison of large cycles and small cycles for different spectral methods.
Method Large cycles Small cycles n l PDF of amplitude n s PDF of amplitude JM ν 0 , P ⋅ T Eq. (22) ν 0 , H ⋅ T Rayleigh FC ν 0 , L ⋅ T Eq. (25) ν 0 , H - ν 0 , L ⋅ T Rayleigh MFC ν 0 , P ⋅ T Eq. (25) ν 0 , H - ν 0 , P ⋅ T RayleighBecause of complexity of the convolution integration, several researches assert that (26) has no analytical solution [9, 10]. To solve this problem, Jiao and Moan used an approximate model (i.e., (22)) to obtain a closed-form solution [8]. However, the approximate model may lead to errors in some cases as in Figure 3 which illustrates the divergence of (19) and (22) for different values of λ0,L∗ and λ0,H∗. It is found that (22) becomes closer to (19) with the increase of λ0,L∗.Figure 3
The divergence of (19) and (22) for JM method in case of (a) λ0,L∗=0.01 and (b) λ0,L∗=0.99.
(a) (b)For FC and MFC methods, (26) was calculated by numerical technique. Although the numerical technique can give a fatigue damage prediction, it is complex and not convenient when applied in real engineering. In addition, the solutions in some cases are not reasonable. In Section 4, an analytical solution of (26) will be derived to evaluate the fatigue damage, and the derivation of the analytical solution focuses on the fatigue damage of the large cycles.


## 3.1. Jiao-Moan (JM) Method

To simplify the study,X(t), XL(t), and XH(t) are normalized as X∗(t), XL∗(t), and XH∗(t) through the following transformation:(14)X∗t=Xtλ0=XLtλ0+XHtλ0=XL∗t+XH∗tand then(15)λ0∗=λ0,L∗+λ0,H∗=1,where(16)λ0,L∗=λ0,Lλ0,λ0,H∗=λ0,Hλ0.Jiao-Moan points out that the small stress cycles are produced by the envelope of the HF process, which follows the Rayleigh distribution. The fatigue damage due to the small stress cycles can be obtained according to (5).While the large stress cycles are from the envelop process,P(t) (see Figure 1), the amplitude of P(t) is equal to(17)Qt=RLt+RHt,where RL(t) and RH(t) are the envelopes of XL∗(t) and XH∗(t), respectively.Figure 1
Bimodal processX∗(t), the envelope process P(t), and the amplitude process Q(t).The distribution ofQ(t) can be written as a form of a convolution integral(18)fQq=∫0qfRLq-xfRHxdx=∫0qfRLyfRHq-ydy.RL(t) and RH(t) obey the Rayleigh distribution; therefore, (18) has an analytical solution which is given [8](19)fQq=qλ0,L∗·exp⁡-q22λ0,L∗+qλ0,H∗·exp⁡-q22λ0,H∗+exp⁡-q22·2πλ0,L∗λ0,H∗·Φqλ0,L∗λ0,H∗+Φqλ0,H∗λ0,L∗-1·q2-1.The rate of mean zero upcrossing due toP(t) can be calculated as(20)ν0,P=λ0,L∗ν0,L1+λ0,H∗λ0,L∗ν0,Hν0,LδH2,where(21)δH=1-λ1,H∗2λ0,H∗λ2,H∗.An approximation was made by Jiao and Moan for (19) as follows [8]:(22)fQq≈λ0,L∗-λ0,L∗λ0,H∗·q·exp⁡-q22λ0,L∗+2πλ0,L∗λ0,H∗·q2-1·exp⁡-q22.After the approximation, a closed-form solution of the bandwidth correction factor can be then derived [8](23)ρ=ν0,Pν0×λ0,L∗m/2+21-λ0,H∗λ0,L∗+πλ0,L∗λ0,H∗mΓm/2+1/2Γm/2+1+ν0,Hν0λ0,H∗m/2.Finally, the fatigue damage can be obtained as (6) and (12).


## 3.2. Fu-Cebon (FC) Method

Similarly to JM method, Fu and Cebon also considered that the total damage is produced by a large cycle (SH+SL) and a small cycle (SL), as depicted in Figure 2. The small cycles are from the HF process, and the distribution of the amplitude PSs(S) is a Rayleigh distribution, as shown in (4). However, the number of cycles associated with the small cycles ns is different from JM method and equals ν0,H-ν0,L·T. According to (5), the damage due to the small cycles is(24)Ds=ν0,H-ν0,L·TK2λ0,HmΓ1+m2.Figure 2
The large cycles and small cycles for a random stress process.The amplitude of the large cyclesSl can be approximated as the sum of amplitude of the LF and HF processes, the distribution of which can be expressed by a convolution of two Rayleigh distributions [9]. (25)PSlS=∫0SPSLyPSHS-ydy=∫0SPSLS-yPSHydy=1λ0,Lλ0,He-S2/2λ0,H∫0SSy-y2e-Uy2+VSydy,where U=1/2λ0,L+1/2λ0,H and V=1/λ0,H.The number of cycles of the large cycles isnl=ν0,L·T. Thus, the fatigue damage due to the large stress cycles can be expressed by(26)Dl=ν0,L·TK∫0∞SmPSlSdS.Equation (26) can be calculated with numerical integration [9, 10]. Therefore, the total damage can be obtained according to (13).


## 3.3. Modify Fu-Cebon (MFC) Method

Benasciutti and Tovo made a comparison between JM method and FC method and concluded that using the envelop process is more suitable [10]. Thus, a hybrid technique is adopted to modify the FC method. More specifically, the large cycles and small cycles are produced according to the idea of FC method. The number of cycles associated with the large cycles is defined similarly to JM method. That is, nl=ν0,P·T, while the number of cycles corresponding to the small cycles is ns=ν0,H-ν0,P·T. The total damage for MFC method can be then written according to (13).Although the accuracy of the MFC method is improved, the fatigue damage still has to be calculated with numerical integral.


## 3.4. Comparison of Three Bimodal Methods

Detailed comparison of the aforementioned three bimodal methods can be found in Table1. In all methods, the amplitude of the small cycle obeys Rayleigh distribution, and the corresponding fatigue damage has an analytical expression as in (5); the distribution of amplitude of the large cycle is convolution integration of two Rayleigh distributions, and the relevant fatigue damage can be calculated by (26).Table 1
Comparison of large cycles and small cycles for different spectral methods.
Method Large cycles Small cycles n l PDF of amplitude n s PDF of amplitude JM ν 0 , P ⋅ T Eq. (22) ν 0 , H ⋅ T Rayleigh FC ν 0 , L ⋅ T Eq. (25) ν 0 , H - ν 0 , L ⋅ T Rayleigh MFC ν 0 , P ⋅ T Eq. (25) ν 0 , H - ν 0 , P ⋅ T RayleighBecause of complexity of the convolution integration, several researches assert that (26) has no analytical solution [9, 10]. To solve this problem, Jiao and Moan used an approximate model (i.e., (22)) to obtain a closed-form solution [8]. However, the approximate model may lead to errors in some cases as in Figure 3 which illustrates the divergence of (19) and (22) for different values of λ0,L∗ and λ0,H∗. It is found that (22) becomes closer to (19) with the increase of λ0,L∗.Figure 3
The divergence of (19) and (22) for JM method in case of (a) λ0,L∗=0.01 and (b) λ0,L∗=0.99.
(a) (b)For FC and MFC methods, (26) was calculated by numerical technique. Although the numerical technique can give a fatigue damage prediction, it is complex and not convenient when applied in real engineering. In addition, the solutions in some cases are not reasonable. In Section 4, an analytical solution of (26) will be derived to evaluate the fatigue damage, and the derivation of the analytical solution focuses on the fatigue damage of the large cycles.


## 4. Derivation of an Analytical Solution

### 4.1. Derivation of an Analytical Solution for (25)

Equation (25) can be rewritten as (27)PSlS=∫0SySλ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy-∫0Sy2λ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy.Equation (27) will be divided into two items.(  1) The First Item. It is as follows:(28)I1=∫0SySλ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy=-Sλ0,L+λ0,Hexp⁡-S22λ0,H+Sλ0,L+λ0,H·exp⁡-S22λ0,L+S2λ0,Lλ0,L+λ0,H·2πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H.(  2) The Second Item. It is as follows:(29)I2=∫0Sy2λ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy=-Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,H-Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,L-2Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,H-exp⁡-S22λ0,L+exp⁡-S22λ0,L+λ0,H2πλ0,Lλ0,Hλ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L1+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H1λ0,L+λ0,H+S2λ0,Hλ0,Lλ0,L+λ0,H2.The analytical solution of (25) can be then obtained(30)PSlS=Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,L+Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,H+S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H.Note that whenλ0,L+λ0,H=1, (30) is just equal to (19) derived by Jiao and Moan [8]. Therefore, (19) is a special case of (30).


### 4.2. Derivation of an Analytical Solution for (26) Based on (30)

The derivation of an analytical solution for (26) is on the basis of (30), as (31)Z=∫0∞Sm·PSlSdS=∫0∞Sm·Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,L+Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,HdS+∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,HdS.Equation (31) will be divided into five parts.(  1) The First Part. It is as follows:(32)Z1=∫0∞Sm·Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,LdS=λ0,L2λ0,L+λ0,H2·2λ0,Lm·Γ1+m2.(  2) The Second Part. It is as follows:(33)Z2=∫0∞Sm·Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,HdS=λ0,H2λ0,L+λ0,H2·2λ0,Hm·Γ1+m2.(  3) The Third Part. It is as follows:(34)Z3=1λ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,H∫0∞Sm+2-Smλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HdS=m2πλ1λ2λ0,L+λ0,Hm-22m-1Γm+12.(  4) The Fourth Part. It is as follows:(35)Z4=∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,HdS.Equation (35) contains a standard Normal cumulative distribution function Φ(·). It is difficult to get an exact solution directly. Thus, a new variable is introduced, as follows:(36)t=Sλ0,L+λ0,H.With a method of variable substitution, (35) can be simplified:(37)Z4=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2∫0∞tm+2exp⁡-t22Φtλ0,Lλ0,Hdt-∫0∞tmexp⁡-t22Φtλ0,Lλ0,Hdt.By defining(38)Hλ0,L,λ0,H,m=∫0∞tmexp⁡-t22Φtλ0,Lλ0,Hdt,(37) becomes(39)Z4=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2Hλ0,L,λ0,H,m+2-Hλ0,L,λ0,H,mand using integration by parts, (38) reduces to(40)Hλ0,L,λ0,H,m=m-1·Hλ0,L,λ0,H,m-2+1212πλ0,Lλ0,H2λ0,Hλ0,L+λ0,HmΓm2.Equation (40) is a recurrence formula; when m is an odd number, it becomes(41)Hλ0,L,λ0,H,m=m-1!!·Hλ0,L,λ0,H,1+λ0,Lλ0,Hm-1!!22π∑k=2m+1/212k-2!!2λ0,Hλ0,L+λ0,H2k-1Γ2k-12,where ·!! is a double factorial function and H(λ0,L,λ0,H,1) has an analytical expression which can be derived conveniently(42)Hλ0,L,λ0,H,1=12+12λ0,Lλ0,L+λ0,H.When m is an even number, H(λ0,L,λ0,H,m) is(43)Hλ0,L,λ0,H,m=m-1!!·Hλ0,L,λ0,H,0+λ0,Lλ0,Hm-1!!22π∑k=1m/212k-1!!2λ0,Hλ0,L+λ0,H2kΓk,where(44)Hλ0,L,λ0,H,0=2π4+12ππ2-atan⁡λ0,Hλ0,L.Specific derivation of (44) can be seen in Appendix A.In addition, a Matlab program has been written to calculateH(λ0,L,λ0,H,m) in Appendix B.(  5) The Fifth Part. It is as follows:(45)Z5=∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Hλ0,Lλ0,L+λ0,HdS.Similarly to the fourth part, the analytical solution of the fifth part can be derived as(46)Z5=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2Hλ0,H,λ0,L,m+2-Hλ0,H,λ0,L,m.The final solution is(47)Z=Z1+Z2-Z3+Z4+Z5.


## 4.1. Derivation of an Analytical Solution for (25)

Equation (25) can be rewritten as (27)PSlS=∫0SySλ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy-∫0Sy2λ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy.Equation (27) will be divided into two items.(  1) The First Item. It is as follows:(28)I1=∫0SySλ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy=-Sλ0,L+λ0,Hexp⁡-S22λ0,H+Sλ0,L+λ0,H·exp⁡-S22λ0,L+S2λ0,Lλ0,L+λ0,H·2πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H.(  2) The Second Item. It is as follows:(29)I2=∫0Sy2λ0,Lλ0,Hexp⁡-y22λ0,H·exp⁡-S-y22λ0,Ldy=-Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,H-Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,L-2Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,H-exp⁡-S22λ0,L+exp⁡-S22λ0,L+λ0,H2πλ0,Lλ0,Hλ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L1+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H1λ0,L+λ0,H+S2λ0,Hλ0,Lλ0,L+λ0,H2.The analytical solution of (25) can be then obtained(30)PSlS=Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,L+Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,H+S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,H.Note that whenλ0,L+λ0,H=1, (30) is just equal to (19) derived by Jiao and Moan [8]. Therefore, (19) is a special case of (30).


## 4.2. Derivation of an Analytical Solution for (26) Based on (30)

The derivation of an analytical solution for (26) is on the basis of (30), as (31)Z=∫0∞Sm·PSlSdS=∫0∞Sm·Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,L+Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,HdS+∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,H-1+ΦSλ0,Hλ0,Lλ0,L+λ0,HdS.Equation (31) will be divided into five parts.(  1) The First Part. It is as follows:(32)Z1=∫0∞Sm·Sλ0,Lλ0,L+λ0,H2exp⁡-S22λ0,LdS=λ0,L2λ0,L+λ0,H2·2λ0,Lm·Γ1+m2.(  2) The Second Part. It is as follows:(33)Z2=∫0∞Sm·Sλ0,Hλ0,L+λ0,H2exp⁡-S22λ0,HdS=λ0,H2λ0,L+λ0,H2·2λ0,Hm·Γ1+m2.(  3) The Third Part. It is as follows:(34)Z3=1λ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,H∫0∞Sm+2-Smλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HdS=m2πλ1λ2λ0,L+λ0,Hm-22m-1Γm+12.(  4) The Fourth Part. It is as follows:(35)Z4=∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Lλ0,Hλ0,L+λ0,HdS.Equation (35) contains a standard Normal cumulative distribution function Φ(·). It is difficult to get an exact solution directly. Thus, a new variable is introduced, as follows:(36)t=Sλ0,L+λ0,H.With a method of variable substitution, (35) can be simplified:(37)Z4=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2∫0∞tm+2exp⁡-t22Φtλ0,Lλ0,Hdt-∫0∞tmexp⁡-t22Φtλ0,Lλ0,Hdt.By defining(38)Hλ0,L,λ0,H,m=∫0∞tmexp⁡-t22Φtλ0,Lλ0,Hdt,(37) becomes(39)Z4=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2Hλ0,L,λ0,H,m+2-Hλ0,L,λ0,H,mand using integration by parts, (38) reduces to(40)Hλ0,L,λ0,H,m=m-1·Hλ0,L,λ0,H,m-2+1212πλ0,Lλ0,H2λ0,Hλ0,L+λ0,HmΓm2.Equation (40) is a recurrence formula; when m is an odd number, it becomes(41)Hλ0,L,λ0,H,m=m-1!!·Hλ0,L,λ0,H,1+λ0,Lλ0,Hm-1!!22π∑k=2m+1/212k-2!!2λ0,Hλ0,L+λ0,H2k-1Γ2k-12,where ·!! is a double factorial function and H(λ0,L,λ0,H,1) has an analytical expression which can be derived conveniently(42)Hλ0,L,λ0,H,1=12+12λ0,Lλ0,L+λ0,H.When m is an even number, H(λ0,L,λ0,H,m) is(43)Hλ0,L,λ0,H,m=m-1!!·Hλ0,L,λ0,H,0+λ0,Lλ0,Hm-1!!22π∑k=1m/212k-1!!2λ0,Hλ0,L+λ0,H2kΓk,where(44)Hλ0,L,λ0,H,0=2π4+12ππ2-atan⁡λ0,Hλ0,L.Specific derivation of (44) can be seen in Appendix A.In addition, a Matlab program has been written to calculateH(λ0,L,λ0,H,m) in Appendix B.(  5) The Fifth Part. It is as follows:(45)Z5=∫0∞Sm·S2-λ0,L+λ0,Hλ0,L+λ0,H22πλ0,Lλ0,Hλ0,L+λ0,Hexp⁡-S22λ0,L+λ0,HΦSλ0,Hλ0,Lλ0,L+λ0,HdS.Similarly to the fourth part, the analytical solution of the fifth part can be derived as(46)Z5=2πλ0,Lλ0,Hλ0,L+λ0,Hm-2Hλ0,H,λ0,L,m+2-Hλ0,H,λ0,L,m.The final solution is(47)Z=Z1+Z2-Z3+Z4+Z5.


## 5. Numerical Validation

In this part, the accuracy of FC method and JM method and the derived analytical solution will be validated with numerical integration. Transformation of (26) will be carried out first.


### 5.1. Treatment of Double Integral Based on FC Method

As pointed out by Fu and Cebon and Benasciutti and Tovo [9, 10], FC’s numerical integration can be calculated as the following processes.Equation (26) contains a double integral, in which variables S and y are in the range of 0,∞ and 0,S, respectively. Apparently, the latter is not compatible with the integration range of Gauss-Legendre quadrature formula. Therefore, by using a integration transformation(48)y=S21+t,the integral part of (26) can be simplified to(49)J=∫0∞SmPSlSdS=1λ0,Lλ0,H∫0∞∫-11SmS231-t2exp⁡-S/221+t22λ0,H·exp⁡-S/221-t22λ0,LdtdS.Equation (49) can be thus calculated with Gauss-Legendre and Gauss-Laguerre quadrature formula.


### 5.2. Treatment of Numerical Integral for (31)

Direct calculation of (31) may lead to some mathematical accumulative errors. To obtain a precise integral solution, (31) has to be handled with a variable substitution, and the result is(50)Z′=λ0,Lm+4/2+λ0,Hm+4/2λ0,L+λ0,H2·2m·Γ1+m2+2πλ0,Lλ0,Hλ0,L+λ0,Hm-2∫0∞tm+2-tmexp⁡-t22Φtλ0,Lλ0,H+Φtλ0,Hλ0,L-1dt.Note that Z1 and Z2 for (31) have analytical solution; therefore, only Z3, Z4, and Z5 are dealt with in (50).The solution of (50) can be obtained through Gauss-Laguerre quadrature formula.The accuracy of the numerical integral in (49) and (50) depends on the order of nodes and weights which can be obtained from handbook of mathematics [25]. The accuracy increases with the increasing orders. However, from the engineering point of view, too many orders of nodes and weights will lead to difficulty in calculation. The integral results of (50) are convergent when the orders of nodes and weights are equal to 30. Therefore, the present study takes the order of 30.


### 5.3. Discussion of Results

It is very convenient to use JM method in the case ofλ0,L+λ0,H=1, while FC method and MFC method can be used for any case regardless of λ0,L and λ0,H. Therefore, the analytical results are divided into two: λ0,L+λ0,H=1 and λ0,L+λ0,H≠1.The solutions calculated by different mathematical methods form=3 and m=4 in the case of λ0,L+λ0,H=1 are plotted in Figures 4 and 5. It turns out that the proposed analytical solution gives the same results with the exact numerical integral solution for any value of λ0,L and λ0,H. FC’s numerical integral solution matches the exact numerical integral solution very well in most cases. However, for relatively low value of λ0,L or when λ0,L tends to 1, it may not be right. JM’s approximate solution is close to the exact numerical integral solution only in a small range.Figure 4
Comparison of different solutions form=3.Figure 5
Comparison of different solutions form=4.FC’s numerical integral solutions, the exact numerical integral solutions, and the analytical solutions form=3 and m=4 in the case of λ0,L+λ0,H≠1 are shown in Tables 2 and 3. The results indicate that the latter two solutions are always approximately equal for any value of λ0,L and λ0,H, while FC’s numerical integral solutions show good agreement with the exact integral solutions only in a few cases. As like the case of λ0,L+λ0,H=1, for relatively high or low values of λ0,L and λ0,H (i.e., λ0,L=100, λ0,H=1×10-5), FC’s integral solutions may give incorrect results. This phenomenon is in accord with the analysis of Benasciutti and Tovo (e.g., FC’s numerical integration may be impossible for too low values of λ0,L and λ0,H) [10].Table 2
Comparison of different solutions form=3 in the case of λ0,L+λ0,H≠1.
λ 0 , L λ 0 , H FC’s numerical integral solution for (49) The exact numerical integral solution for (50) The analytical solution for (47) 1.000 E - 05 9.000 E - 05 3.657 E - 12 6.183 E - 06 6.183 E - 06 3.000 E - 03 7.000 E - 03 3.292 E - 03 7.591 E - 03 7.591 E - 03 5.000 E - 02 5.000 E - 02 2.591 E - 01 2.522 E - 01 2.522 E - 01 2.000 E - 01 6.000 E - 01 5.262 E + 00 5.267 E + 00 5.267 E + 00 5.000 E + 00 4.000 E + 00 2.146 E + 02 2.146 E + 02 2.146 E + 02 8.000 E + 01 2.000 E + 01 7.062 E + 03 7.062 E + 03 7.062 E + 03 4.000 E + 02 6.000 E + 02 1.144 E + 05 2.493 E + 05 2.493 E + 05 5.000 E + 03 3.000 E + 03 1.238 E + 04 5.602 E + 06 5.602 E + 06 6.000 E + 07 5.000 E + 07 8.468 E - 05 9.180 E + 12 9.180 E + 12 7.000 E + 09 4.000 E + 09 9.073 E - 09 8.999 E + 15 8.999 E + 15 1.000 E + 02 1.000 E - 05 3.084 E - 01 3.762 E + 03 3.762 E + 03Table 3
Comparison of different solutions form=4 in the case of λ0,L+λ0,H≠1.
λ 0 , L λ 0 , H FC’s numerical integral solution for (49) The exact numerical integral solution for (50) The analytical solution for (47) 1.000 E - 05 9.000 E - 05 2.580 E - 13 1.437 E - 07 1.437 E - 07 3.000 E - 03 7.000 E - 03 1.187 E - 03 1.832 E - 03 1.832 E - 03 5.000 E - 02 5.000 E - 02 2.198 E - 01 1.942 E - 01 1.942 E - 01 2.000 E - 01 6.000 E - 01 1.131 E + 01 1.130 E + 01 1.130 E + 01 5.000 E + 00 4.000 E + 00 1.567 E + 03 1.567 E + 03 1.567 E + 03 8.000 E + 01 2.000 E + 01 1.682 E + 05 1.682 E + 05 1.682 E + 05 4.000 E + 02 6.000 E + 02 6.632 E + 06 1.915 E + 07 1.915 E + 07 5.000 E + 03 3.000 E + 03 7.706 E + 05 1.216 E + 09 1.216 E + 09 6.000 E + 07 5.000 E + 07 5.308 E - 03 2.344 E + 17 2.344 E + 17 7.000 E + 09 4.000 E + 09 5.688 E - 07 2.289 E + 21 2.289 E + 21 1.000 E + 02 1.000 E - 05 7.230 E - 01 8.006 E + 04 8.006 E + 04In a word, for any value ofλ0,L and λ0,H, the analytical solution derived in this paper always gives an accurate result. Furthermore, it can be solved very conveniently and quickly with the aid of a personal computer through a program given in Appendix B.


## 5.1. Treatment of Double Integral Based on FC Method

As pointed out by Fu and Cebon and Benasciutti and Tovo [9, 10], FC’s numerical integration can be calculated as the following processes.Equation (26) contains a double integral, in which variables S and y are in the range of 0,∞ and 0,S, respectively. Apparently, the latter is not compatible with the integration range of Gauss-Legendre quadrature formula. Therefore, by using a integration transformation(48)y=S21+t,the integral part of (26) can be simplified to(49)J=∫0∞SmPSlSdS=1λ0,Lλ0,H∫0∞∫-11SmS231-t2exp⁡-S/221+t22λ0,H·exp⁡-S/221-t22λ0,LdtdS.Equation (49) can be thus calculated with Gauss-Legendre and Gauss-Laguerre quadrature formula.


## 5.2. Treatment of Numerical Integral for (31)

Direct calculation of (31) may lead to some mathematical accumulative errors. To obtain a precise integral solution, (31) has to be handled with a variable substitution, and the result is(50)Z′=λ0,Lm+4/2+λ0,Hm+4/2λ0,L+λ0,H2·2m·Γ1+m2+2πλ0,Lλ0,Hλ0,L+λ0,Hm-2∫0∞tm+2-tmexp⁡-t22Φtλ0,Lλ0,H+Φtλ0,Hλ0,L-1dt.Note that Z1 and Z2 for (31) have analytical solution; therefore, only Z3, Z4, and Z5 are dealt with in (50).The solution of (50) can be obtained through Gauss-Laguerre quadrature formula.The accuracy of the numerical integral in (49) and (50) depends on the order of nodes and weights which can be obtained from handbook of mathematics [25]. The accuracy increases with the increasing orders. However, from the engineering point of view, too many orders of nodes and weights will lead to difficulty in calculation. The integral results of (50) are convergent when the orders of nodes and weights are equal to 30. Therefore, the present study takes the order of 30.


## 5.3. Discussion of Results

It is very convenient to use JM method in the case ofλ0,L+λ0,H=1, while FC method and MFC method can be used for any case regardless of λ0,L and λ0,H. Therefore, the analytical results are divided into two: λ0,L+λ0,H=1 and λ0,L+λ0,H≠1.The solutions calculated by different mathematical methods form=3 and m=4 in the case of λ0,L+λ0,H=1 are plotted in Figures 4 and 5. It turns out that the proposed analytical solution gives the same results with the exact numerical integral solution for any value of λ0,L and λ0,H. FC’s numerical integral solution matches the exact numerical integral solution very well in most cases. However, for relatively low value of λ0,L or when λ0,L tends to 1, it may not be right. JM’s approximate solution is close to the exact numerical integral solution only in a small range.Figure 4
Comparison of different solutions form=3.Figure 5
Comparison of different solutions form=4.FC’s numerical integral solutions, the exact numerical integral solutions, and the analytical solutions form=3 and m=4 in the case of λ0,L+λ0,H≠1 are shown in Tables 2 and 3. The results indicate that the latter two solutions are always approximately equal for any value of λ0,L and λ0,H, while FC’s numerical integral solutions show good agreement with the exact integral solutions only in a few cases. As like the case of λ0,L+λ0,H=1, for relatively high or low values of λ0,L and λ0,H (i.e., λ0,L=100, λ0,H=1×10-5), FC’s integral solutions may give incorrect results. This phenomenon is in accord with the analysis of Benasciutti and Tovo (e.g., FC’s numerical integration may be impossible for too low values of λ0,L and λ0,H) [10].Table 2
Comparison of different solutions form=3 in the case of λ0,L+λ0,H≠1.
λ 0 , L λ 0 , H FC’s numerical integral solution for (49) The exact numerical integral solution for (50) The analytical solution for (47) 1.000 E - 05 9.000 E - 05 3.657 E - 12 6.183 E - 06 6.183 E - 06 3.000 E - 03 7.000 E - 03 3.292 E - 03 7.591 E - 03 7.591 E - 03 5.000 E - 02 5.000 E - 02 2.591 E - 01 2.522 E - 01 2.522 E - 01 2.000 E - 01 6.000 E - 01 5.262 E + 00 5.267 E + 00 5.267 E + 00 5.000 E + 00 4.000 E + 00 2.146 E + 02 2.146 E + 02 2.146 E + 02 8.000 E + 01 2.000 E + 01 7.062 E + 03 7.062 E + 03 7.062 E + 03 4.000 E + 02 6.000 E + 02 1.144 E + 05 2.493 E + 05 2.493 E + 05 5.000 E + 03 3.000 E + 03 1.238 E + 04 5.602 E + 06 5.602 E + 06 6.000 E + 07 5.000 E + 07 8.468 E - 05 9.180 E + 12 9.180 E + 12 7.000 E + 09 4.000 E + 09 9.073 E - 09 8.999 E + 15 8.999 E + 15 1.000 E + 02 1.000 E - 05 3.084 E - 01 3.762 E + 03 3.762 E + 03Table 3
Comparison of different solutions form=4 in the case of λ0,L+λ0,H≠1.
λ 0 , L λ 0 , H FC’s numerical integral solution for (49) The exact numerical integral solution for (50) The analytical solution for (47) 1.000 E - 05 9.000 E - 05 2.580 E - 13 1.437 E - 07 1.437 E - 07 3.000 E - 03 7.000 E - 03 1.187 E - 03 1.832 E - 03 1.832 E - 03 5.000 E - 02 5.000 E - 02 2.198 E - 01 1.942 E - 01 1.942 E - 01 2.000 E - 01 6.000 E - 01 1.131 E + 01 1.130 E + 01 1.130 E + 01 5.000 E + 00 4.000 E + 00 1.567 E + 03 1.567 E + 03 1.567 E + 03 8.000 E + 01 2.000 E + 01 1.682 E + 05 1.682 E + 05 1.682 E + 05 4.000 E + 02 6.000 E + 02 6.632 E + 06 1.915 E + 07 1.915 E + 07 5.000 E + 03 3.000 E + 03 7.706 E + 05 1.216 E + 09 1.216 E + 09 6.000 E + 07 5.000 E + 07 5.308 E - 03 2.344 E + 17 2.344 E + 17 7.000 E + 09 4.000 E + 09 5.688 E - 07 2.289 E + 21 2.289 E + 21 1.000 E + 02 1.000 E - 05 7.230 E - 01 8.006 E + 04 8.006 E + 04In a word, for any value ofλ0,L and λ0,H, the analytical solution derived in this paper always gives an accurate result. Furthermore, it can be solved very conveniently and quickly with the aid of a personal computer through a program given in Appendix B.


## 6. Case Study

In this section, the bandwidth correction factor is used to compare different bimodal spectral methods.A general analytical solution (GAS) of fatigue damage for JM, FC, and MFC method can be written as(51)DGAS=nlK·Z+nsK2λ0,HmΓ1+m2.Z can be obtained according to (47), and nl and ns can be chosen as defined in Table 1, which represent different bimodal spectra methods, that is, JM, FC, and MFC method.According to (6), the analytical solution of the bandwidth correction factor is(52)ρGAS=DGASDNB.Likewise, a general integration solution (GIS) as given in FC and MFC method can be written as(53)DGIS=nlK·J+nsK2λ0,HmΓ1+m2.J can be obtained from (49).The integration solution of the damage correction factor is(54)ρGIS=DGISDNB.


### 6.1. Ideal Bimodal Spectra

Different spectral shapes have been investigated in the literatures [5, 26]. Bimodal PSDs with two rectangular blocks are used to carry out numerical simulations, as shown in Figure 6. Two blocks are characterized by the amplitude levels SωL and SωH, as well as the frequency ranges ωb-ωa and ωd-ωc.Figure 6
Ideal bimodal spectra.ω L and ωH are the central frequencies, as defined in(55)ωL=ωa+ωb2,ωH=ωc+ωd2.A 1 and A2 represent the areas of block 1 and block 2 and are equal to the zero-order moment, respectively. For convenience, the sum of the two spectral moments is normalized to unity; that is, A1+A2=1.To ensure that these two spectra are approximately narrow band processes,ωb/ωa=ωd/ωc=1.1.Herein, two new parameters are introduced(56)B=λ0,Hλ0,L,R=ωHωL.The parameter values in the numerical simulations are conducted as follows:B=0.1, 0.4, 1, 2, and 9; R=6 and 10 which ensure two spectra are well-separated; m=3, 4, and 5 which are widely used in welded steel structures. ωa does not affect the simulated bandwidth correction factor ρ [10, 27]. Thus, ωa can be taken as the arbitrary value in the theory as long as ωa>0; herein, ωa=5rad/s.In the process of time domain simulation, the time series generated by IDFT [28] contains 20 million sample points and 150–450 thousand rainflow cycles, which is a sufficiently long time history, so that the sampling errors can be neglected.Figure7 is the result of JM method. The bandwidth correction factor calculated by (23) is in good agreement with the results obtained from (52). For m=3, JM method can provide a reasonable damage prediction. However, for m=4 and m=5, this method tends to underestimate the fatigue damage.Figure 7
The bandwidth correction factor of JM method for (a)R=6 and (b) R=10.
(a) (b)Figures8 and 9 display the results of FC method and MFC method, which are both calculated by the proposed analytical solution (see (52)) and FC’s numerical integration solution (see (54)). The bandwidth correction factor calculated by (52) is very close to the results obtained from (54). Besides, the FC method always provides conservative results compared with RFC. The MFC method improves the accuracy of the original FC method to some extent. However, it may underestimate the fatigue damage in some cases.Figure 8
The bandwidth correction factor of FC method for (a)R=6 and (b) R=10.
(a) (b)Figure 9
The bandwidth correction factor of MFC method for (a)R=6 and (b) R=10.
(a) (b)As has been discussed above, the bandwidth correction factor calculated by the analytical solution has divergence with that computed from RFC. The disagreement is not because of the error from the analytical solution but arises from the fact that the physical models of the original bimodal methods and the rainflow counting method are different.


### 6.2. Real Bimodal Spectra

In practice, the rectangular bimodal spectra in previous simulations cannot represent real spectra encountered in the structures. Therefore, more realistic bimodal stress spectra will be chosen to predict the vibration-fatigue-life by using the analytical solution proposed in this paper. For offshore structures subjected to random wave loading, the PSDs of fatigue stress of joints always exhibit two predominant peaks of frequency. Wirsching gave a general expression to characterize the PSD, and the model has been widely used in a few surveys [10, 23, 29]. The analytical form is(57)Sω=AHsexp⁡-1050/TD4ω4TD4ω41-ω/ωn22+2ζω/ωn2,where A is a scale factor, Hs is the significant wave height, TD is the dominant wave period, ωn is the natural angular frequency of structure, and ζ is the damping coefficient.A and Hs do not affect the value of the bandwidth correction factor ρ [10, 27]. Thus, they are chosen to be equal to unity for simplicity. The value of other parameters can be seen in Table 4.Table 4
Parameters for real bimodal spectra.
A H s (m) T D (s) ω N (rad/s) ζ Group 1 1 1 15 2.2 1.5% Group 2 1 1 5 10 1%Real stress spectra corresponding to group 1 and group 2 can be seen in Figure10. S(ω) is a double peak spectrum, the first peak is produced by the peak of random wave spectrum, and the second one is excitated by the first mode response of structures.Figure 10
Real bimodal spectra corresponding to (a) group 1 and (b) group 2.
(a) (b)The results of different bimodal methods are shown in Figures11(a) and 11(b). The bandwidth correction factor for JM method is calculated with (23) and (52). It should be noted that λ0,L+λ0,H≠1 for real bimodal spectra in Figures 10(a) and 10(b), so, λ0,L and λ0,H should be normalized as (14), (15), and (16), while ρ for FC method and MFC method is obtained through (52) and (54).Figure 11
The bandwidth correction factor of different bimodal methods for (a) group 1 and (b) group 2.
(a) (b)Under the real stress spectra, the zero-order spectral moments corresponding to low frequency and high frequency are very small. JM method provides acceptable damage estimation compared with RFC. The results computed by (54) for FC and MFC method may be not correct in some cases. The incorrect results are mainly caused by (49). However, the proposed analytical solution (see (52)) can give a satisfactory damage prediction and always provide a conservative prediction. More importantly, (52) is more convenient to apply in predicting vibration-fatigue-life than numerical integration.


## 6.1. Ideal Bimodal Spectra

Different spectral shapes have been investigated in the literatures [5, 26]. Bimodal PSDs with two rectangular blocks are used to carry out numerical simulations, as shown in Figure 6. Two blocks are characterized by the amplitude levels SωL and SωH, as well as the frequency ranges ωb-ωa and ωd-ωc.Figure 6
Ideal bimodal spectra.ω L and ωH are the central frequencies, as defined in(55)ωL=ωa+ωb2,ωH=ωc+ωd2.A 1 and A2 represent the areas of block 1 and block 2 and are equal to the zero-order moment, respectively. For convenience, the sum of the two spectral moments is normalized to unity; that is, A1+A2=1.To ensure that these two spectra are approximately narrow band processes,ωb/ωa=ωd/ωc=1.1.Herein, two new parameters are introduced(56)B=λ0,Hλ0,L,R=ωHωL.The parameter values in the numerical simulations are conducted as follows:B=0.1, 0.4, 1, 2, and 9; R=6 and 10 which ensure two spectra are well-separated; m=3, 4, and 5 which are widely used in welded steel structures. ωa does not affect the simulated bandwidth correction factor ρ [10, 27]. Thus, ωa can be taken as the arbitrary value in the theory as long as ωa>0; herein, ωa=5rad/s.In the process of time domain simulation, the time series generated by IDFT [28] contains 20 million sample points and 150–450 thousand rainflow cycles, which is a sufficiently long time history, so that the sampling errors can be neglected.Figure7 is the result of JM method. The bandwidth correction factor calculated by (23) is in good agreement with the results obtained from (52). For m=3, JM method can provide a reasonable damage prediction. However, for m=4 and m=5, this method tends to underestimate the fatigue damage.Figure 7
The bandwidth correction factor of JM method for (a)R=6 and (b) R=10.
(a) (b)Figures8 and 9 display the results of FC method and MFC method, which are both calculated by the proposed analytical solution (see (52)) and FC’s numerical integration solution (see (54)). The bandwidth correction factor calculated by (52) is very close to the results obtained from (54). Besides, the FC method always provides conservative results compared with RFC. The MFC method improves the accuracy of the original FC method to some extent. However, it may underestimate the fatigue damage in some cases.Figure 8
The bandwidth correction factor of FC method for (a)R=6 and (b) R=10.
(a) (b)Figure 9
The bandwidth correction factor of MFC method for (a)R=6 and (b) R=10.
(a) (b)As has been discussed above, the bandwidth correction factor calculated by the analytical solution has divergence with that computed from RFC. The disagreement is not because of the error from the analytical solution but arises from the fact that the physical models of the original bimodal methods and the rainflow counting method are different.


## 6.2. Real Bimodal Spectra

In practice, the rectangular bimodal spectra in previous simulations cannot represent real spectra encountered in the structures. Therefore, more realistic bimodal stress spectra will be chosen to predict the vibration-fatigue-life by using the analytical solution proposed in this paper. For offshore structures subjected to random wave loading, the PSDs of fatigue stress of joints always exhibit two predominant peaks of frequency. Wirsching gave a general expression to characterize the PSD, and the model has been widely used in a few surveys [10, 23, 29]. The analytical form is(57)Sω=AHsexp⁡-1050/TD4ω4TD4ω41-ω/ωn22+2ζω/ωn2,where A is a scale factor, Hs is the significant wave height, TD is the dominant wave period, ωn is the natural angular frequency of structure, and ζ is the damping coefficient.A and Hs do not affect the value of the bandwidth correction factor ρ [10, 27]. Thus, they are chosen to be equal to unity for simplicity. The value of other parameters can be seen in Table 4.Table 4
Parameters for real bimodal spectra.
A H s (m) T D (s) ω N (rad/s) ζ Group 1 1 1 15 2.2 1.5% Group 2 1 1 5 10 1%Real stress spectra corresponding to group 1 and group 2 can be seen in Figure10. S(ω) is a double peak spectrum, the first peak is produced by the peak of random wave spectrum, and the second one is excitated by the first mode response of structures.Figure 10
Real bimodal spectra corresponding to (a) group 1 and (b) group 2.
(a) (b)The results of different bimodal methods are shown in Figures11(a) and 11(b). The bandwidth correction factor for JM method is calculated with (23) and (52). It should be noted that λ0,L+λ0,H≠1 for real bimodal spectra in Figures 10(a) and 10(b), so, λ0,L and λ0,H should be normalized as (14), (15), and (16), while ρ for FC method and MFC method is obtained through (52) and (54).Figure 11
The bandwidth correction factor of different bimodal methods for (a) group 1 and (b) group 2.
(a) (b)Under the real stress spectra, the zero-order spectral moments corresponding to low frequency and high frequency are very small. JM method provides acceptable damage estimation compared with RFC. The results computed by (54) for FC and MFC method may be not correct in some cases. The incorrect results are mainly caused by (49). However, the proposed analytical solution (see (52)) can give a satisfactory damage prediction and always provide a conservative prediction. More importantly, (52) is more convenient to apply in predicting vibration-fatigue-life than numerical integration.


## 7. Conclusion

In this paper, bimodal spectral methods to predict vibration-fatigue-life are investigated. The conclusions are as follows:(1)
An analytical solution of the convolution integral of two Rayleigh distributions is developed. Besides, this solution is a general form which is different from that proposed by Jiao and Moan. The latter is only a particular case.(2)
An analytical solution based on bimodal spectral methods is derived. It is validated that the analytical solution shows a good agreement with numerical integration. More importantly, the analytical solution has a stronger attraction than numerical integration in engineering application.(3)
For JM method, the original approximate solution (see (23)) is reasonable in most cases; the analytical solution (see (52)) can give better prediction.(4)
In ideal bimodal processes, JM method and MFC method may overestimate or underestimate the fatigue damage, while in real bimodal processes, they give a conservative prediction. FC method always provides conservative results in any cases. Therefore, FC method can be recommended as a safe design technique in real engineering structures.



---

*Source: 1010726-2017-01-31.xml*