Life estimation of the beam with normal distribution parameters and subjected to cyclic load

1. Introduction

The beam is the common structural element of the civil and mechanical engineering. It could withstand the bending moment and shear force primarily. If the beam failure occurs during service, it should lead to the huge economic lost or the casualties. The failure due to fatigue might be the common one. In mechanical engineering, there are a lot of factors which result in the beam with random parameters. The factors should include manufacturing error of machine tool, ambient temperature, random vibration of the foundation of machine tool, material uniformity and so on. Moreover, many beams used in mechanical engineering commonly carry the cyclic load, such as the loads due to all kinds of faults of the rotating machine. Therefore, it is important to estimate the life of the beam with random parameters and subjected to the cyclic load.

The random parameters of the beam result in the random stress under the deterministic cyclic load. The known methods of estimating the life of metallic structures and components under random stress can be classified into two groups: time-domain method based on a cycle counting method and frequency-domain method using power spectral density (PSD) of stochastic processes [1]. One of time-domain method, the combination of rainflow counting and Palmgren-Miner rule, is recognized as the “gold standard” for estimating fatigue damage from variable amplitude stress or for comparison to other estimate techniques [2, 3]. However, it has been proved to spend a lot of time in processing data when time-domain method is employed to estimate the life under random stress [4, 5]. For the reason, frequency-domain method has been developed. It can produce (exact or approximated) analytical expressions by which we can estimate cycle distribution and fatigue life under a given counting procedure, without, however, requiring knowledge of each stress or strain time history [6, 7]. The life estimation of the beam under random stress has reported by some documents. The focus of the reported works is to predict the life of the beam with the deterministic parameters and withstanding the random loading. For example, the fatigue life prediction method of the dual cantilever beam subjected to the narrow band random load was introduced by Rajnikant Bhogilal Thakkar [8] and was based on the fatigue life under the variable sinusoidal load. The fatigue life of the Aluminum 6061 T6 cantilever beam, where the input was the white noise acceleration, were predicted using seven types of the frequency methods and Rainflow counting method respectively and they were compared with these obtained by experiments in [9]. The effect of the damping ratios of the test structure on the fatigue life was also analyzed by the four types of tests with the different damping ratios. An analytic methodology was developed by Paulus et al. [10] to predict the fatigue life of the cantilever beam experiencing stationary and Gaussian random vibration excitation. In addition, Paulus et al. [11] proposed a semi-empirical life method of the cantilevered beam subjected to the random excitation with a wide range of PSD (power spectral density) shapes and amplitudes. Two methods were employed by Lin et al. [12, 13] to predict the fatigue lives of RC beams strengthened with CFL (carbon fiber laminate) subjected to the exciting force due to the random traffic flow. They were based on the residual life model and the residual strength model respectively. Moreover, the fatigue life of the beam with the uncertain parameters was discussed by a few works. For example, the effect of the uncertainties on the fatigue life of a higher order sandwich beam was discussed by Wang and Sun [14]. In their work, the uncertainties included the density of foam core materials and the thicknesses of the core and the top bottom skins, two independent white noise loadings with the same spectral density function were applied to two different points on the top skin.

It has not been reported to estimate the fatigue life of the beam which has several random parameters and is subjected to cyclic load. Therefore, it should be discussed in this work. Here, these random parameters are assumed to be normally distributed. The remainder of this paper is organized as follows: In Section 2, it should be tested by Monte-Carlo simulation that the generalized displacement and velocity follow normal distributions with the periodic mean and variance estimated by the random perturbation method, and the mean, variance and correlation coefficient of the random dynamic stress and its time derivative are formulated using the first-order approximation in a Taylor series. In Section 3, the expected cumulative damage based on Palmgren-Miner rule is modeled where the random dynamic stress is non-stationary and follows the normal distribution with the nonzero mean at any time $t$. The conclusions are drawn finally.

2. Random dynamical stress analysis of the beam with normal distribution parameters

2.1. Random dynamical response analysis

When the beam undergoes small random vibrations around a position of stable equilibrium and is assumed to be with viscous damping, the dynamical equation of the beam with random parameters is:

1

$\mathbf{M}\ddot{\mathbf{X}}+\mathbf{C}\dot{\mathbf{X}}+\mathbf{K}\mathbf{X}=\mathbf{F}\left(t\right),$

where $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are the $n$×$n$ dimension mass matrix, damping matrix and stiffness matrix respectively and $\mathbf{F}\left(t\right)$ is the $n$ dimension external load vector of which the element is the cyclic external load. $\ddot{\mathbf{X}}$, $\dot{\mathbf{X}}$ and $\mathbf{X}$ are the $n$ dimension generalized acceleration, velocity and displacement vector respectively. $n$ is the total number of the degree of freedom. One or multiple elements of $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ are the random variables. All random parameters of the considered beam are written by one random vector $\mathbf{B}$. Generally, $\mathbf{B}$ is composed of the total length $L$, width $b$, height $h$, modulus of elasticity $E$ and density $\rho$ which are supposed to be normally distributed and mutually independent. The random perturbation method [15, 16] is employed to analyze the random dynamical response of the dynamical Eq. (1). The mean of the generalized displacement vector is:

2

$E\left[\mathbf{X}\left(t\right)\right]=E\left[{\mathbf{X}}_d\left(t\right)\right]={\mathbf{X}}_d\left(t\right),$

where ${\mathbf{X}}_d\left(t\right)$ is the deterministic part of the generalized displacement vector $\mathbf{X}$ and is solved by:

3

${\mathbf{M}}_d{\ddot{\mathbf{X}}}_d+{\mathbf{C}}_d{\dot{\mathbf{X}}}_d+{\mathbf{K}}_d{\mathbf{X}}_d=\mathbf{F}\left(t\right),$

where ${\mathbf{M}}_d$, ${\mathbf{C}}_d$ and ${\mathbf{K}}_d$ are the deterministic part of $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$ respectively. They are obtained by substituting the mean of all random parameters into $\mathbf{M}$, $\mathbf{C}$ and $\mathbf{K}$. The covariance matrix of the generalized displacement vector is:

4

$\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathbf{X}\right)=\left[\begin{array}{cccc}\mathrm{V}\mathrm{a}\mathrm{r}\left(x\left(1\right)\right)& & & \mathrm{s}\mathrm{y}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{y}\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(x\left(2\right),x\left(1\right)\right)& \mathrm{V}\mathrm{a}\mathrm{r}\left(x\left(2\right)\right)& & \\ ⋮& ⋮& \ddots & \\ \mathrm{C}\mathrm{o}\mathrm{v}\left(x\left(n\right),x\left(1\right)\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(x\left(n\right),x\left(2\right)\right)& \cdots & \mathrm{V}\mathrm{a}\mathrm{r}\left(x\left(n\right)\right)\end{array}\right]={\mathbf{X}}_{\mathbf{B}}\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathbf{B}\right){\mathbf{X}}_{\mathbf{B}}^{\mathrm{T}},$

where ${\mathbf{X}}_{\mathbf{B}}$ is the sensitivity matrix of ${\mathbf{X}}_d\left(t\right)$ to the random vector $\mathbf{B}$ and is:

5

${\mathbf{X}}_{\mathbf{B}}=\left[\begin{array}{lll}{x}_{d,1}\left(1\right)& x_{d,2}\left(1\right)\cdots & x_{d,m}\left(1\right)\\ \begin{array}{l}{x}_{d,1}\left(2\right)\\ ⋮\end{array}& \begin{array}{l}{x}_{d,2}\left(2\right)\cdots \\ ⋮\end{array}& \begin{array}{l}{x}_{d,m}\left(2\right)\\ ⋮\end{array}\\ x_{d,1}\left(n\right)& x_{d,2}\left(n\right)\cdots & x_{d,m}\left(n\right)\end{array}\right].$

$\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathbf{B}\right)$ is the covariance matrix of the random parameter vector $\mathbf{B}$ and is:

6

$\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathbf{B}\right)=\left[\begin{array}{cccc}\mathrm{V}\mathrm{a}\mathrm{r}\left(b_1\right)& & & \mathrm{s}\mathrm{y}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{y}\\ \mathrm{C}\mathrm{o}\mathrm{v}\left(b_2,b_1\right)& \mathrm{V}\mathrm{a}\mathrm{r}\left(b_2\right)& & \\ ⋮& ⋮& \ddots & \\ \mathrm{C}\mathrm{o}\mathrm{v}\left(b_m,b_1\right)& \mathrm{C}\mathrm{o}\mathrm{v}\left(b_m,b_2\right)& \cdots & \mathrm{V}\mathrm{a}\mathrm{r}\left(b_m\right)\end{array}\right],$

where $m$ is the total number of random parameters and $x_{d,j}\left(i\right)=\partial x_d\left(i\right)/\partial b_j$ which is solved by:

7

${\mathbf{M}}_d{\ddot{\mathbf{X}}}_{d,j}+{\mathbf{C}}_d{\dot{\mathbf{X}}}_{d,j}+{\mathbf{K}}_d{\mathbf{X}}_{d,j}=-\left({\mathbf{M}}_{d,j}{\ddot{\mathbf{X}}}_d+{\mathbf{C}}_{d,j}{\dot{\mathbf{X}}}_d+{\mathbf{K}}_{d,j}{\mathbf{X}}_d\right).$

${\mathbf{M}}_{d,\mathrm{}\mathrm{}j}={\left.\frac{\partial \mathbf{M}}{\partial b_j}\right|}_{\mathbf{B}={\mathbf{B}}_d},{\mathbf{C}}_{d,\mathrm{}\mathrm{}j}={\left.\frac{\partial \mathbf{C}}{\partial b_j}\right|}_{\mathbf{B}={\mathbf{B}}_d},{\mathbf{K}}_{d,\mathrm{}\mathrm{}j}={\left.\frac{\partial \mathbf{K}}{\partial b_j}\right|}_{\mathbf{B}={\mathbf{B}}_d},$

where ${\mathbf{B}}_d$ is the deterministic part of $\mathbf{B}$.

The skewness and kurtosis of the generalized displacement and velocity of node 1, 3 and 5 of the random cantilever beam at 0.005 and 0.012 second are shown in table 1. All random parameters’ $CV$ is the identical number in one Monte-Carlo simulation. The samples of the generalized displacement are generated by Monte-Carlo method and Newmark-$\beta$ method. The sample size is 5000. The considered cantilever beam with normal distribution parameters is shown in Fig. 1(b). It is divided into 5 elements. The means of the total length $L$, width $b$, height $h$, density ${\rho }_1$ and elastic modulus $E$ are 0.25 m, 0.02 m, 0.03 m, 7850 kg/m³ and 2.06×10¹¹ Pa respectively. The deterministic cyclic load $F\left(t\right)=3500\mathrm{s}\mathrm{i}\mathrm{n}\left(100\pi t\right)$ is applied to the free end node. From Table 1, it could be trusted that the random generalized displacement and velocity at any time are the normal distribution when all CVs of the random parameters are small (generally $CV\le$0.01).

The displacement mean and standard deviation of node 1 of the cantilever beam with normal distribution parameters calculated by Monte-Carlo simulation and random perturbation method respectively with $CV=$0.01 are shown in Fig. 2 and 3. Here, the values of the parameters except $CV$ are identical with the numerical calculation of the skewness and kurtosis of the generalized displacement and velocity. The difference between Monte-Carlo simulation and random perturbation method is very little according to Fig. 2 and 3. In fact, the similar results are found in the mean and standard deviation of the generalized displacement and velocity of other nodes. Therefore, it might be one solution to estimating the mean and variance by random perturbation method when the random generalized displacement and velocity of the beam with normal distribution parameters are normally distributed at any time.

Table 1Skewness and kurtosis of generalized displacement and velocity of node 1, 3 and 5 of cantilever beam with normal distribution random parameters at 0.005 and 0.012 second

$CV$	Time (s)		Node 1				Node 3				Node 5
			$x\left(1\right)$	$\dot{x}\left(1\right)$	$x\left(2\right)$	$\dot{x}\left(2\right)$	$x\left(5\right)$	$\dot{x}\left(5\right)$	$x\left(6\right)$	$\dot{x}\left(6\right)$	$x\left(9\right)$	$\dot{x}\left(9\right)$	$x\left(10\right)$	$\dot{x}\left(10\right)$
0.005	0.005	Skewness	0.060	0.072	0.057	0.067	0.060	0.071	0.057	0.066	0.060	0.071	0.057	0.066
		Kurtosis	3.007	3.010	3.004	3.007	3.007	3.009	3.004	3.007	3.007	3.009	3.004	3.007
	0.012	Skewness	0.063	–0.053	0.059	0.052	0.062	–0.054	0.059	–0.052	0.062	–0.054	0.059	–0.053
		Kurtosis	3.007	3.006	3.005	3.003	3.007	3.006	3.005	3.003	3.007	3.006	3.005	3.003
0.010	0.005	Skewness	0.148	0.171	0.145	0.168	0.148	0.170	0.145	0.166	0.148	0.169	0.145	0.165
		Kurtosis	3.019	3.027	3.073	3.070	3.019	3.027	3.073	3.070	3.019	3.026	3.073	3.070
	0.012	Skewness	0.153	–0.135	0.150	–0.131	0.153	–0.136	0.150	-0.132	0.153	–0.137	0.149	–0.134
		Kurtosis	3.021	3.017	3.072	3.084	3.020	3.017	3.073	3.082	3.020	3.017	3.073	3.081
0.050	0.005	Skewness	0.744	0.882	0.712	0.834	0.744	0.876	0.713	0.825	0.744	0.870	0.713	0.817
		Kurtosis	4.138	4.541	4.056	4.413	4.134	4.522	4.057	4.383	4.135	4.563	4.058	4.362
	0.012	Skewness	0.776	–0.681	0.740	–0.662	0.775	–0.686	0.738	-0.668	0.773	–0.690	0.737	–0.673
		Kurtosis	4.229	3.999	4.140	3.936	4.225	4.009	4.133	3.952	4.220	4.019	4.128	3.963

Fig. 1a) Coordinate system and b) finite element model of the beam

a)

b)

Fig. 2Displacement mean of node 1 of the considered cantilever beam with CV= 0.01

Fig. 3Displacement standard deviation of node 1 of the considered cantilever beam with CV= 0.01

3. Life estimation of a beam with normal distribution parameters and subjected to cyclic load

According to Palmgren-Miner rule, the expected cumulative damage during the time interval $t_1$ to $t_2$ is [19]:

where $\mathrm{C}$ and $m$ are constants depending on the material, ${\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}(y,\mathbf{B},t)$ and ${\dot{\sigma }}_{z\mathrm{m}\mathrm{a}\mathrm{x}}(y,\mathbf{B},t)$ are replaced by ${\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\dot{\sigma }}_{z\mathrm{m}\mathrm{a}\mathrm{x}}$ respectively for simplicity, ${\stackrel{-}{v}}_{{\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}}^{+}\left(t\right)$ is the expected number of crossings (at level ${\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}$, with positive slope) and is calculated by [20]:

If ${\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}$ and ${\dot{\sigma }}_{z\mathrm{m}\mathrm{a}\mathrm{x}}$ are normally distributed, their joint probability density function is:

where $\rho =\rho (y,t)$ is estimated by Eq. (19), ${\stackrel{-}{\sigma }}_{z\mathrm{m}\mathrm{a}\mathrm{x}}={\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}(y,{\mathbf{B}}_d,t)$ and ${\stackrel{-}{\dot{\sigma }}}_{z\mathrm{m}\mathrm{a}\mathrm{x}}={\dot{\sigma }}_{z\mathrm{m}\mathrm{a}\mathrm{x}}(y,{\mathbf{B}}_d,t)$. Then, the expected number of crossings is:

where $\mathrm{E}\mathrm{r}\mathrm{f}\mathrm{c}\left(z\right)=1-\mathrm{e}\mathrm{r}\mathrm{f}\left(z\right)$, erf(z) is the error function defined by:

and:

Fig. 4Mean of the maximal fixed-end stress of the considered cantilever beam when CV= 0.005 and other parameters are the same as Section 2.1

Fig. 5Standard deviation of the maximal fixed-end stress of the considered cantilever beam when CV= 0.005 and other parameters are the same as Section 2.1

Fig. 6Correlation coefficient between maximal fixed-end stress and its time derivative of the considered cantilever beam when CV= 0.005 and other parameters are the same as Section 2.1

Then, the partial derivative of the expected number of crossings with respect to ${\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}$ is:

where:

Then, the expected cumulative damage of the beam with normal distribution parameters and subjected to cyclic load during the time interval t₁ to t₂ could be calculated by substituting Eq. (25) and (26) into Eq. (20) and using the numerical integration method.

Using the proposed method, the expected instantaneous damages (${\int }_0^{\mathrm{\infty }}\partial {\stackrel{-}{v}}_{{\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}}^{+}\left(t\right)/\partial {\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}{\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}^{m}d{\sigma }_{z\mathrm{m}\mathrm{a}\mathrm{x}}$) of the cantilever beam with normal distribution parameters are shown in Fig. 7 and 8 when $CV=$0.005 and 0.010 respectively. Other parameters of the considered beam are the same as Section 2.1. The expected instantaneous damage of the discussed random cantilever beam subjected to cyclic load is the periodic function of the time $t$ and its amplitude with $CV=$0.010 is more than one with $CV=$0.005 from Fig. 7 and 8. The expected cumulative damage with $CV=$0.005 and 0.010 during one second are 4.945×10^-4 and 1.742×10^-3 and their corresponding lives (numbers of cycles to failure) are 1.072×10⁵ and 3.043×10⁴ respectively.

Fig. 7Expected instantaneous damage of the cantilever beam with normal distribution parameters and subjected to cyclic load when CV= 0.005 and other parameters are the same as Section 2.1

Fig. 8Expected instantaneous damage of the cantilever beam with normal distribution parameters and subjected to cyclic load when CV= 0.01 and other parameters are the same as Section 2.1

4. Conclusions

One method for estimating the life of the beam with normal distribution parameters and subjected to cyclic load was presented in this work. The random generalized displacement and velocity of the beam with random parameters at any time are the normal distributions when CVs of all random parameters are small (generally $CV\le$0.01). Their means and variances could be accurately estimated by the random perturbation method. By compared with Monte-Carlo simulation, the means, variances and their correlation coefficient of the random dynamic stress and its time derivative could be estimated using the first-order approximation method in a Taylor series. The presented method could estimate the life of the beam with normal distribution parameters and subjected to cyclic load where the random dynamic stress is non-stationary and follows the normal distribution with the nonzero mean at any time $t$.