Stochastic analysis of the critical velocity of an axially moving cracked elastic plate

In this study, a probabilistic analysis of the critical velocity for an axially moving cracked elastic and isotropic plate is presented. Axially moving materials are commonly used in modelling of manufacturing processes, like paper making and plastic forming. In such systems, the most serious threats to runnability are instability and material fracture, and ﬁnding the critical value of velocity is essential for eﬃciency. In this paper, a formula for the critical velocity is derived under constraints for the probabilities of instability and fracture. The signiﬁcance of randomness in diﬀerent model parameters is investigated for parameter ranges typical of paper material and paper machines. The results suggest that the most signiﬁcant factors are variation in the crack length and tension magnitude.


Introduction
In industry, there are many systems the behaviour of which can be described by the mathematical model of an axially moving material. Thus, during the last few decades, the mechanics of such materials have aroused much interest among researchers. Traditionally, the studies of axially moving materials are based on a deterministic approach, although in reality, the problem parameters are not known deterministically. In industrial paper manufacturing, which is one of the application areas of axially moving materials, uncertainty factors include, e.g., the strength of the paper web, variation of tension with respect to space and time in the press system, and defects, which vary in their geometry and location in the web. These factors are considerable: according to Uesaka [1], the majority of web breaks in paper production are caused by tension variations, combined with strength variations of the paper web. Wathén [2] discusses the effect of flaws of paper on web breaks and, according to him, even a seemingly flawless paper can fail at very low tensions due to stress concentrations caused by discontinuities, e.g. cuts and shives, in structure.
Finding the optimal value of velocity for an axially moving material is essential, when the efficiency of the corresponding manufacturing process is considered. The most critical threats to good runnability of such a system are instability and material fracture, and on these phenomena a change in tension magnitude has opposite effects. An increase in tension has a stabilizing effect [3], but high tension may lead to growing or arising of cracks. Web tension too low or too high may cause a web break, which deteriorates production efficiency.

2
The modelling of vibrations of travelling elastic materials has interested many researchers. The first paper on the subject dates from 1897, when Skutch published a paper [4] concerning the axially moving string. The first papers in English were published in the 1950's, when Sack [5] and Archibald and Emslie [6] studied the axially moving string model. Since then, many researchers have continued the studies of moving elastic material. E.g., Wickert and Mote [7] studied stability of axially moving strings and beams using modal analysis and Greens function method. The stability of travelling twodimensional rectangular membranes and plates has been studied, e.g., by Lin [8] and Banichuk et al. [3]. A more extensive literature review of the history of the studies concerning deterministic elastic models can be found in [3], the results of which we also exploit in this study. In the recent studies concerning axially moving plates, material properties such as orthotropicity [9,10] or viscoelasticity [11,12] have been taken into consideration and their effects on the plate behavior have been investigated.
In addition, there are studies considering stationary plates with random parameters. For example, the free transverse vibrations of elastic rectangular plates with random material properties were considered and statistical characteristics of the random eigenvalues were determined by Sobczyk [13]. Wood and Zaman [14] considered a collection of elastic rectangular plates with random inhomogeneities vibrating freely under simply supported boundary conditions. Soares [15] considered uncertainty modelling of plates subjected to compressive loads.
The field of fracture mechanics was developed by Irwin [16], based on the early papers of Inglis [17], Griffith [18] and Westergaard [19]. Various deterministic analysis of vibrations and stability of stationary cracked beams and plates exist in the literature. For a literature review, we refer to [20].
As far as the authors know, axially moving materials have been studied in a stochastic setup only in [21,22]. In these studies, the critical velocity of an axially moving plate was derived in the case in which there is a random length crack on the plate, or the tension, to which the plate is subjected, varies randomly. This research extends these studies by introducing several other parameters as random variables simultaneously in the model. In this paper, we also compare the effect of introducing variation between different problem parameters, in order to decide the randomness of which parameters is the most significant in terms of the critical velocity. For the analysis, we have chosen the setup and parameter ranges to be applicable for paper material and paper making.
The formula for the critical velocity of the plate is derived under constraints for instability and fracture. Depending on the distributions of the problem parameters, numerical methods may be needed in solving the critical velocity. In the paper industrial example, simultaneously introducing several problem parameters as random leads to the use of numerical methods. Due to its simplicity and accuracy, we use Monte Carlo simulation to solve the problem with several random variables.

Critical velocity of a travelling plate
We consider a rectangular part of an elastic and isotropic band, which is moving at a constant velocity V 0 between supporting rollers. Denoting the part as where and b are prescribed parameters of length and width, the plate is assumed to travel in the x direction. The supporting rollers are located at x = 0 and x = . (See Figure 1.) The considered part D is represented as a thin elastic plate having constant thickness h, Poisson ratio ν, Young modulus E, and bending rigidity The mass of the plate per unit area is denoted by m. It is further assumed that the plate is subjected to homogeneous tension T acting in the x direction.
We consider the case in which there is a single crack in the plate. The length of the crack is denoted by ξ. (See also Figure 1.)

Characterization of instability of the plate
We first briefly present a deterministic stability analysis for a travelling plate without a crack. Especially, we are interested in critical regimes, where the plate approaches its maximum stable velocity. Details of the analysis can be found in [3].
We perform a standard dynamic analysis (see, e.g., [23]). The transverse displacement of the travelling plate is described by the deflection function w, which depends on the space coordinates x, y and time t. It is assumed that the absolute values of the deflection function w and its derivatives are small.
The Kirchhoff plate theory is applied. To study the dynamic behavior of the plate, the following equation for the travelling plate is used: where C = T m and (4) As boundary conditions, the classical simply-supported and free boundary conditions [24,25] are used. The simply supported boundary conditions read as (w) x=0, = 0, and the equations for the boundaries free of tractions can be presented as follows: The solution of the dynamic boundary value problem of (3) -(8) can be represented as whereω is the frequency of small transverse vibrations ands = iω is the complex characteristic parameter;s = Res + i Ims.
If the parameters is purely imaginary andω is real, the plate performs harmonic vibrations of a small amplitude and its motion can be considered stable. If the real part ofs becomes positive, the transverse vibrations grow exponentially and, consequently, the behaviour is unstable (See Figure 2).
It can be shown by dynamic analysis that the travelling plate undergoes divergence instability at a sufficiently high speed [8] and thus it is sufficient 6 to perform static analysis, i.e., study the case withs = 0. (See also Figure   2.) The stationary equations for W are, substituting (9) into (3) and setting with boundary conditions (6) - (8). We rewrite (10) as where is the eigenvalue.
In order to determine the minimal eigenvalue (12) of the problem (6) - (11), and the corresponding eigenfunction W = W (x, y), the following representation is applied: where f (y/b) is an unknown function. Note that W in (13) satisfies the boundary condition (6). The solution, a half-sine in the longitudinal direction, is well-known (see, e.g., [8]). Using the dimensionless formulation and the relations (7) - (8) and (10) - (14), the following eigenvalue problem for the unknown function f (η) is obtained: The problem (15) - (17) is now considered as a spectral boundary value problem. The symmetric solution of (15) - (17) corresponding to the minimal eigenvalue λ is where A and B are arbitrary constants.
In [3] it was shown that the buckling mode is symmetric. When γ ≤ 1, the symmetric divergence mode is thus To determine the constants A and B in (18), we insert f from equation (18) into the boundary conditions (16) - (17): 8 The condition for a non-trivial solution to exist in the form (18) - (19) is that the determinant of the system (20) -(21) must vanish. This leads to the transcendental equation which determines the eigenvalues λ = γ 2 as an implicit function. In (22), and Detailed analysis of the properties of Φ and Ψ can be found in [10].
It follows that the critical velocity of instability of the travelling plate is given by where γ * is the root of the equation (22). The motion of the plate is stable when its velocity satisfies As the root γ * does not depend on the value of tension T , the limit of stable velocity (25) is increased when increasing the tension T . However, in the case of a cracked plate, increasing tension may lead to the growth of the crack.

Probabilistic analysis of critical velocity considering instability and fracture
In the following, we include uncertainty in the parameters of the plate and formulate the problem of finding the critical velocity of the plate un-9 der constraints for the probabilities of instability and fracture. Introducing a constraint for the probability of failure is a way to formulate statistical mechanical problems, used by many researchers (see, e.g., [26,27]).
Let (Ω, F , P) be a probability space, where Ω is a sample space, F is its σ-algebra and P is a probability measure on F . Let ξ, h, m, ν, E, the critical strain energy release rate, G C , and the solution of the equation (22), γ * , be random variables on Ω, the image of each random variable being an interval of {x ∈ R : x > 0}. The tension applied to the supported edges of the plate is defined as where T 0 is a positive constant and θ is a random variable on Ω with The stress intensity factors K j , j = I, II, related to the crack modes I (opening) and II (in-plane shear) have the form and are random variables. In (29), α j is a weight function, the expression of which depends on the crack geometry. Formulae for weight functions α j are presented, e.g., in [28], [29], [30] and [31]. In this study, the weight functions α j are assumed to be strictly positive, continuous and increasing.
The fracture toughness K C of the material is given by and is also a random variable.
Considering both instability and fracture, the problem of critical velocity of the plate in the presence of a mode I crack reads as where p i , p f denote the upper bounds of the probabilities of instability and fracture, respectively. In the case of a mixed mode crack, the constraint (32) is replaced by For the failure criterion (34) of a mixed mode crack, see [32,33].
In the case of a mode I crack we denote The inequalities (32) and (34) are equivalent to where F X is the cumulative distribution function of the random variable X, 11 The function F X is continuous and strictly increasing, when and Thus the maximal value of T 0 that satisfies (32) (or (34)), denoted by T cr 0 , is found on the interval in (40) by solving the equation and is the p f th order quantile of F X : The inequality (33) is equivalent to where The cumulative distribution function F Y T 0 is continuous, increasing with respect to V 0 and decreasing with respect to T 0 . Hence the maximum value of V 0 , such that (33) and (32) (or (34)) are satisfied, is found by studying The solution of the problem (31)-(33) (or the problem (31), (33) and (34)) , denoted by V cr 0 , is the p i th order quantile of F Y : Table 1 shows the solutions of the problems in which only the crack length, the tension variation and the strain energy release rate or one of them is regarded as a random variable. When the only random valued parameter is ξ or G C , there are no random variables in the constraint (33). In this case, the critical velocity is simply given by (25). In Table 1, is the cumulative distribution function of the random variable j.

Numerical approximation of the quantiles
and the pth quantile F −1 (p) can be approximated as where np is the first integer ≥ np. In [34] it is shown that S (n) np is a weakly consistent quantile estimator.
With the sample S 1 , . . . , S n from the distribution F , the strongly consis- where χ A is the indicator of the event A, and its standard error is The estimator (50) and its standard error (51) are used for evaluating the sufficient sample size in quantile estimation. By setting where r ∈ (0, 1), a lower bound for the sample size n in estimating the pth The values of the quantiles F −1 θ and F −1 G C , with normally distributed θ and G C , can be obtained by a subroutine of a standard program library (e.g. SciPy). For a Weibull distributed ξ, an analytical expression for F −1 ξ exists.

Numerical results and discussion
In this section, the solutions obtained in Sections 2.2 and 2.3 are illustrated in a paper industrial context.
The random variables ξ, h, ν, E, G C and θ were assumed to be independent. The random variables h, ν, E, G C , θ, were assumed to obey the normal distribution with the means µ j , j = h, ν, E, G C , θ, in Table 2. The mean values were chosen to correspond to the properties of printing paper. The mass of the plate m was assumed to be dependent of h: 14 The effect of variation in the parameter values on the solution was studied by introducing the coefficient of variation, CV j , for the random variable j as a relation of its standard deviation, σ j , and mean, µ j : The crack length was assumed to be distributed according to the Weibull distribution with a shape parameter s ∈ (0, 1]. The cumulative distribution function of the Weibull distributed crack length is where x ≥ 0 and the scale parameter c > 0. The expected value and variance of the crack length are where Γ is the gamma function, When the shape parameter s ∈ (0, 1], the probability density function is strictly decreasing. For s ∈ (0, 1), the probability density function tends to infinity when x tends to zero. For s = 1, the probability density function tends to 1/c when x tends to zero, the distribution corresponding to the exponential distribution with intensity 1/c. Figure 3 shows the effect of changing the value of c on the probability density function. It shows that when c is decreased, the probability mass becomes more concentrated near the origin, corresponding to the small values of the crack length. When c is decreased, both the expected value and variance of the crack length decrease. The shape parameter s has the opposite effect on the probability density function, as is seen in Figure 3.
In cases where the problem parameters m, h, ν, E and G C were regarded as constants, their values were set to the mean values in Table 2. If the crack length was regarded as a constant, its value was set to be equal to the is not a significant factor, when the critical velocity of the cracked thin plate is studied.
In Figure 5 the critical velocities given by the model with random G C , θ and ξ and the models with a single random parameter are compared to the critical velocity given by the corresponding fully deterministic model. The critical velocities were thus divided by the coefficientṼ 0 , γ * is the solution of the equation (22) with ν = µ ν and κ = /πµ b , The crack length distribution parameters were set as s = 0.5 and c = 0.005.
The coefficientṼ 0 is the critical velocity of the fully deterministic model with a crack of the length E[ξ] = 0.01 (m). The relative critical velocity of the fully deterministic model becomes equal to unity. Figure 5 shows that, with a fixed reliability level, among the models with a single random valued parameter, the model with random crack length produces the lowest critical velocity. The model with random strain energy release rate provides critical values that are the closest to the critical values by the deterministic model. It is seen that, with low CV G C , variation in strain energy release rate may not be significant, the size of the effect depending on the probability p f . The effect of variation in tension is more significant than the impact of variation in strain energy release rate.
The lowest critical velocity in Figure 5 is given by the model with random strain energy release rate, tension and crack length. The effect of CV θ and CV G C on the relative critical velocity was studied in the range CV θ , CV G C ∈ [0, 0.2], with fixed crack length distribution (s = 0.5, c = 0.005). A change in CV G C had no effect on the relative critical velocity but increasing CV θ decreases the critical velocity, although this effect may not be significant with low CV θ .
The above results suggest that, when the effect of variation in the problem parameters is considered, the most significant factor is the random crack length. Including also variation of tension in the model may lower the critical velocity notably, at least if the variation is remarkable. The results agree with the study by Uesaka [1] to some extent. According to [1], the majority of web breaks in paper production are caused by tension variations, combined with strength variations of the paper web. The results of the present study suggest that variation in tension is significant, but variation in strain energy release rate, which causes variation in paper fracture toughness, is not.   Table 3. The obtained critical velocities are higher than the running speeds of the current paper machines. However, it should be noted that the surrounding air is excluded from the model. The presence of surrounding air is known to influence the critical velocity [35,36]. According to [35], the critical velocity obtained with the vacuum model may be even four times the critical value by the model that includes the surrounding air.
The parameter values were chosen to correspond to those of a dry paper web. If wet paper is modelled, also material properties such as viscoelasticity should be included in the model. However, the effect of viscoelasticity on the critical stable velocity is found small (see [11,12]). Another excluded material property, independent of the moisture of paper, is orthotropicity. As with viscoelasticity, the effect of orthotropicity on the critical stable velocity is found small (see [3,10,37]). From the viewpoint of the application, another notable simplification considers the tension profile. In this study, the profile of tension was assumed to be homogeneous, although in paper machines, the measured tension varies in the cross direction [38].

Conclusions
In this study, stochastic analysis of the critical velocity of an axially moving elastic and isotropic plate in the presence of a crack was presented.
The critical velocity was derived under constraints for the probabilities of instability and fracture, which are the most critical threats to runnability of axially moving materials.
The solution was illustrated with a paper industrial example. The crack length was assumed to obey the Weibull distribution while the other random valued problem parameters were assumed to obey the normal distribution.
Altering the set of random parameters in the model, the effect of variation in problem parameters on the critical velocity was studied. The critical 20 velocity obtained by the model where the random parameters included the mass and thickness of the plate, the elastic coefficients, the strain energy release rate, the tension to which the plate is subjected, and the crack length, was compared with the model with only random strain energy release rate, tension and crack length. The results showed no significant effect of introducing variation in the values of the mass and thickness of the plate or in the elastic coefficients.
When the models in which only the strain energy release rate, tension or crack length is regarded as a random variable were compared, the lowest The results suggest that when the critical velocity of the cracked plate is studied, it is essential to include variation in crack length, and possibly also variation in tension magnitude, in the model. Randomness of other problem parameters is not significant, and the model can be simplified such that the random parameters include only the crack length and tension.