Probabilistic analysis of a cantilever beam subjected to random loads via probability density functions

Abstract

This paper addresses the probabilistic analysis of the deflection of a cantilever beam by means of a randomization of the classical governing fourth-order differential equation with null boundary conditions. The probabilistic study is based on the calculation of the first probability density function of the solution, which is a stochastic process, as well as the density function of further quantities of interest associated with this engineering problem such as the maximum slope and deflection at the free end of the cantilever beam, that are treated as random variables. In addition, the probability density function of the bending moment and the shear force will also be computed. The study takes extensive advantage of the so called Random Variable Transformation method, also known as Probability Transformation Method, that allows us to fully unify the probabilistic analysis in three relevant cases commonly studied in the deterministic setting. All the theoretical findings are illustrated via detailed numerical examples corresponding to each one of the three scenarios.

1 Introduction

Deterministic data models based on differential equations have been widely applied in Engineering to design and analyse from simple to complex civil engineering structures. This classical approach guarantees certain structural safety satisfying standard prescriptions. The deterministic analysis involves quantifying the model parameters, which, in practice, are determined after experimentation. The collected data required to set the values of model parameters often contain errors coming from measuring devices (epistemic uncertainty). Apart from this source of uncertainty, the heterogeneity of materials, the lack of knowledge of complex physical phenomena affecting the surrounding medium of the structure including the unknown flow of vehicles, people, etc., are also factors that should be considered to properly calibrate the model parameters (aleatoric uncertainty). These reasons motivate that over the last decades numerous probabilistic-based methods had been proposed to better design and analyse civil engineering structures. The probabilistic approach has crucial advantages such as quantifying structural safety and reliability that cannot be assessed using the deterministic methods. Probabilistic methods avoid the use of over-simplified rules of thumb and permit better understanding the real risk of engineering structures.

The aim of this paper is to contribute to the development and application of stochastic methods to study foundational problems in the realm of civil engineering structures. We shall specifically deal with the stochastic analysis of the deflection (i.e., the amount of displacement, or sag) experienced by a load-carrying beam. In our subsequent analysis we will consider a horizontal beam on which a vertical force is acting. All loaded beams will deflect to a greater or lesser degree, depending upon factors such as the size and placement of loads, the manner of supporting the beam, the beam material and the stiffness of the beam. These two latter factors are rarely known in a deterministic manner because of the heterogeneity of materials, so stochastic methods are more realistic to better analyse the deflection of load-carrying beam.

In the case of analysis of the static deflection of a beam, it is well-known that this phenomenon can be mathematically described by means of the following fourth-order differential equation (Öchsner 2021, Entry 3 of Table 2.4)

$$\begin{aligned} \frac{{\textrm{d}}^4 Y(x)}{{\textrm{d}} x^4}=\frac{W(x)}{Ei}, \quad 0<x<l, \end{aligned}$$

(1)

where Y(x) represents the deflection curve of the beam, E is the Young’s modulus of elasticity of the material of the beam, i denotes the moment of inertia of the cross section of the beam around a horizontal line through the centroid of the cross section, l is the length of the beam and W(x) represents the density of downward force acting vertically on the beam at the spatial point x,  which can be interpreted as load per unit length acting on the beam.

Despite model (1) is usually applied assuming a deterministic (nominal) value for the Young’s modulus, E,  it depends on the physical material properties with which the beam has been built. Due to the heterogeneity of the material, the mathematical nature of the Young’s modulus is random rather than deterministic as it has been reported in different investigations (Khakurel et al. 2021; Soyarslan et al. 2019; Hien et al. 2020; Farsi et al. 2017). On the other hand, the nature of W(x) is also stochastic since, as it has been previously indicated, represents the density of downward force acting perpendicularly on the beam. At every spatial point x,  the value of W(x) will depend on the heterogeneity of the material carried on the beam. Additionally, in practice the value of W(x) on the whole spatial domain is approximated via measurements, so containing epistemic errors. All these considerations lead to rigorously treat (1) as a random differential equation, where E is a random variable and W(x) is a stochastic process defined in a common complete probability space $(\Omega , {\mathcal {F}}, {\mathbb {P}})$. As it shall be seen later, we will assume that E is an absolutely continuous random variable, so having a probability density function (p.d.f.), while for W(x) we will consider both the case that is a parametric stochastic process, which depends on absolutely continuous random variables, and when it is a non-parametric process.

We point out that throughout the paper, and following the standard notation in Probability Theory, we will distinguish random variables and their realizations (which are deterministic quantities), by using in the notation capital and lowercase letters, respectively.

Assuming that the beam is fixed at the origin, i.e., $Y(0)=Y'(0)=0,$ and adding adequate endpoint values, model (1) allows us to study the deflection of beams with different supports: simply supported ($Y(l)=Y''(l)=0$), built-in or fixed end ($Y(l)=Y'(l)=0$) and free end ($Y''(l)=Y'''(l)=0$), for example. In this paper we deal with the stochastic analysis of this latter case, which corresponds to a cantilever (a beam firmly fastened at the origin, $x=0,$ but with no support at the end, $x=l$). The stochastic analysis of simply supported or fixed end beams can also be done using the same probabilistic technique that we shall present throughout the paper. Therefore, hereinafter the following boundary conditions will be assumed

$$\begin{aligned} Y(0)=0,\quad Y'(0)=0,\quad Y''(l)=0,\quad Y'''(l)=0. \end{aligned}$$

(2)

Under our approach, the solution of the randomized boundary problem (1)–(2) is a stochastic process, Y(x),  and our goal will be to determine, under very general assumptions on the model parameters, E and W(x),  the so-called first p.d.f. (1-p.d.f.) of the solution, $f_{Y(x)}(y)$ (Soong 1973, Ch. 3). The calculation of this function is of great usefulness from a practical standpoint since it permits determining the mean, $\mu _Y(x),$ and the variance, $\sigma _Y^2(x),$ of the deflection at each spatial point x of the beam,

$$\begin{aligned} \mu _Y(x)={\mathbb {E}}[Y(x)] = \int _{-\infty }^{\infty } y f_{Y(x)}(y)\,{\textrm{d}}y, \end{aligned}$$

(3)

and

$$\begin{aligned} \sigma _Y^2(x)={\mathbb {E}}[(Y(x)-\mu _Y(x))^2]={\mathbb {E}}[(Y(x))^2]- (\mu _Y(x))^2, \end{aligned}$$

(4)

respectively. Here ${\mathbb {E}}[\cdot ]$ denotes the expectation operator. Additionally, further one-dimensional higher moments with respect to (w.r.t.) the origin and centred at the mean can also be calculated by means of the 1-p.d.f.,

$$\begin{aligned} {\mathbb {E}}[(Y(x))^m] = \int _{-\infty }^{\infty } y^m f_{Y(x)}(y)\,{\textrm{d}}y,\quad m=1,2,\ldots \end{aligned}$$

and

$$\begin{aligned} {\mathbb {E}}[(Y(x)-\mu _Y(x))^m]=\int _{-\infty }^{\infty } (y-\mu _Y(x))^m f_{Y(x)}(y)\,{\textrm{d}}y, \quad m=1,2,\ldots , \end{aligned}$$

respectively. Fixed an arbitrary spatial point x,  the 1-p.d.f. is also particularly relevant to calculate the probability that the deflection, Y(x),  lies on a certain interval of specific interest, say $[y_1,y_2],$

$$\begin{aligned} {\mathbb {P}}[y_1 \le Y(x) \le y_2] = \int _{y_1}^{y_2} f_{Y(x)}(y) {\textrm{d}}y,\quad 0<x < l. \end{aligned}$$

In the stochastic analysis of civil engineering structures this information is of paramount importance since it permits evaluating the probability of breaking of a structure (like a bridge, a balcony, a crane, etc.) that is carrying on distributed or uniform forces. Indeed, this is done by calculating the probability that the deflection of a beam lies within a certain safety interval.

Throughout this paper, we will obtain an explicit expression of the 1-p.d.f., $f_{Y(x)}(y),$ of the random boundary problem (1)–(2) by considering different forms of the density downward force, W(x). To tool to conduct our probabilistic analysis is the Random Variable Transformation (RVT) method, also called Probability Transformation Method (PTM) which allows us to determine the p.d.f. of a random vector that results from transforming another random vector whose p.d.f. is known.

Theorem 1

(RVT method) (Soong 1973, page 25) Let us consider ${\textbf{U}}=(U_1,\ldots ,U_n)$ and ${\textbf{V}}=(V_1,\ldots ,V_n)$ two n-dimensional continuous random vectors defined on a complete probability space $(\Omega ,{\mathcal {F}},{\mathbb {P}}).$ Let ${\textbf{r}}: {\mathbb {R}}^n \rightarrow {\mathbb {R}}^n$ be a one-to-one deterministic transformation of ${\textbf{U}}$ into ${\textbf{V}},$ i.e.,  ${\textbf{V}}={\textbf{r}}({\textbf{U}})$. Assume that ${\textbf{r}}$ is continuous in ${\textbf{U}}$ and has continuous partial derivatives w.r.t. each $U_i,$ $1\le i \le n$. Then,  if $f_{{\textbf{U}}}({\textbf{u}})$ denotes the joint p.d.f. of random vector ${\textbf{U}},$ and ${\textbf{s}}={\textbf{r}}^{-1}=(s_1(v_1,\ldots ,v_n),\ldots ,s_n(v_1,\ldots ,v_n))$ represents the inverse mapping of ${\textbf{r}}=(r_1(u_1,\ldots ,u_n),\ldots ,r_n(u_1,\ldots ,u_n)),$ the joint p.d.f. of random vector ${\textbf{V}}$ is given by

$$\begin{aligned} f_{{\textbf{V}}}({\textbf{v}})=f_{{\textbf{U}}}\left( {\textbf{s}}({\textbf{v}})\right) \left| {\mathcal {J}} \right| , \end{aligned}$$

where $\left| {\mathcal {J}} \right| ,$ which is assumed to be different from zero,  is the absolute value of the Jacobian,  that is defined by the determinant

$$\begin{aligned} {\mathcal {J}}=\det \left( \frac{\partial {\textbf{s}}}{\partial {\textbf{v}}}\right) = \det \left( \begin{array}{ccc} \frac{\partial s_1(v_1,\ldots ,v_n)}{\partial v_1} &{} \cdots &{} \frac{\partial s_n(v_1,\ldots , v_n)}{\partial v_1}\\ \vdots &{} \ddots &{} \vdots \\ \frac{\partial s_1(v_1,\ldots , v_n)}{\partial v_n} &{} \cdots &{} \frac{\partial s_n(v_1,\ldots , v_n)}{\partial v_n}\\ \end{array} \right) . \end{aligned}$$

As it shall be seen later, to apply this key probabilistic result, first we will take advantage of the Laplace transform to explicitly obtain a solution of model (1)–(2) for three different forms of W(x),  hereinafter identified with cases I–III, that are often considered in the analysis of beams in civil engineering.

The analysis of a cantilever beam has been extensively faced, mainly in the deterministic context (see for example, Rizos et al. 1990; Zhang et al. 2022; Le and Ta 2021). In the stochastic setting, the number of contributions is still limited with regard to the specific study of a cantilever beam, although many other types of engineering structures have been studied taking into account uncertainties. Next, we concentrate on commenting the three main approaches, namely, polynomial chaos expansions, Monte Carlo simulation and stochastic finite elements method, that have been applied to analyse a cantilever beam subject to randomness. The next discussion is presented in connection with the study performed throughout this paper. In Pryse and Adhikari (2021), one combines polynomial chaos with Neumann expansion method to obtain closed-form expressions for the first two response moments in different engineering problems. Authors apply the theoretical results obtained to analyze a cantilever beam, where the bending rigidity of the beam, EI,  is assumed to be a stationary Gaussian random field with an exponential-type autocovariance kernel. In Reppel et al. (2019), it is assumed that all parameters of the classical Euler–Bernoulli cantilever beam are random variables. The study is performed in the case of a single load with constant (deterministic) value W. In that paper one also studies the case that the Young’s modulus is a discrete random field modelled first, via simple random walk with stationary independent increments, and, secondly, via an adapted autoregressive model with normal distributed random variables. In both cases, authors use Monte Carlo simulations to numerically compute the distribution of the maximum deflection of a cantilever beam. In Korzeniowski and Weinberg (2019), authors compare stochastic and data-driven finite element methods to study a cantilever by assuming that the Young’s modulus is modelled via a the fluctuation of a nominal (mean) value using a standard Gaussian random variable. In Blondeel et al. (2018), multilevel Monte Carlo is combined with a finite element solver to compute the statistical quantities of the static deflection and frequency response function for a cantilever beam with uncertainty in Young’s modulus under a static and a dynamic load. Authors show that Multilevel Monte Carlo method provides a significant computational cost reduction in comparison with the standard (crude) Monte Carlo method. We also mention the use of Bayesian techniques to better account for the deflection of a cantilever in the case that the spatially variable flexibility is described via a random field represented by means of a Karhunen–Loève expansion (Uribe et al. 2017; Lord et al. 2014). It is worthy to mention that in some contributions, the cantilever beams subject to randomness factors have not been directly analysed, but as specific examples to test new stochastic techniques. In this spirit, we here point out (Wu et al. 2019), where authors propose a stochastic dynamic load identification algorithm to analyse uncertain dynamic system with correlated random system parameters. The adopted method is based on the approximation the Green’s function by the second-order perturbation method and orthogonal polynomial chaos bases. Authors then apply the theoretical results to perform numerical simulations and experimental studies with a cantilever beam under a concentrate stochastic force to estimate the statistical characteristics of the stochastic load from the stochastic structural response samples. To the best of our knowledge, there is still important lack of information with regard the complete stochastic analysis of a cantilever beam subject to uncertainties. This extended study should include the computation of the 1-p.d.f. of the deflection in the case that the most important parameters of the Euler–Bernoulli model (1) are random variables with arbitrary probability distributions and also when the density of the downward force W(x) may randomly vary according to different realistic situations. This paper is aimed at face these aspects from a very general standpoint including the computation of the p.d.f. of relevant quantities associated to the analysis of the cantilever beam as the maximum deflection and the maximum slope at the free end of the beam.

The paper is organized as follows. In Sect. 2, we deal with the case I, corresponding to the scenario that W(x) represents two uniform loads characterized by two independent random variables covering the same span on the beam. In Sect. 3, we study the case II, where multiple independent concentrated loads along the span of the beam are analysed. In this case, W(x) will be represented by means of a train of Dirac delta functions whose intensities (loads) are independent random variables. In Sect. 4, we study the third scenario, case III, where W(x) models a distributed load with an uncertain value characterized by a constant (nominal) value affected by fluctuating changes due to the heterogeneity of the material at every spacial point. Spatial fluctuations are modelled via a stochastic process. In these three above-mentioned sections, we will extensively apply the RVT technique to obtain an explicit expression of the 1-p.d.f. of the solution stochastic process of model (1)–(2). The theoretical results obtained in Sects. 2–4, will be fully illustrated with examples in Sect. 5. In the examples, we will also obtain the p.d.f. of other quantities of interest, such as the maximum deflection and the maximum slope at the free end of the beam, since these quantities are of major interest when analysing and designing cantilever beams in civil engineering satisfying specific preset safety rules. To complete the information of interest in engineering applications, we will also compute the p.d.f.’s of the bending moment and the shear force. Conclusions will be drawn in Sect. 6.

Finally, we point out that the results in the three aforementioned cases I–III are going to be intentionally presented in an analogous manner for the ease of reading, but also to show the generality of the RVT method despite each scenario is completely different from a mathematical standpoint.

2 Case I: Two loads modelled via random variables

This section is addressed to compute the 1-p.d.f., $f_{Y(x)}(y),$ of the stochastic solution of the random model (1)–(2), where

$$\begin{aligned} W(x)= \left\{ \begin{array}{ll} W_0, &{}\quad 0<x \le l/2,\\ W_1, &{}\quad l/2<x \le l. \end{array} \right. \end{aligned}$$

(5)

This model represents the deflection of a cantilever beam with two different loads, represented by $W_0$ and $W_1,$ that occupy the same space on the beam, which is fully loaded. As discussed in the previous section, due to uncertainties associated to data, we will assume that model parameters E,  $W_0$ and $W_1$ are random variables. Hereinafter, $f_{E,W_0,W_1}(e,w_0,w_1)$ will denote the joint p.d.f. of the random vector $(E,W_0,W_1),$ that can be factorized as the product of the corresponding marginal p.d.f., i.e. $f_E(e)f_{W_0}(w_0)f_{W_1}(w_1),$ in the particular, but important case that E,  $W_0$ and $W_1$ are mutually independent. In Fig. 1 we show a graphical scheme of the model.

Fig. 1

Case I: Cantilever beam with two loads modelled via random variables

At this point, it is interesting to note that in our analysis arbitrary distributions have been assumed for model parameters which will provide more generality to our study.

To obtain the 1-p.d.f., $f_{Y(x)}(y),$ taking advantage of the RVT technique (see Theorem 1), we first need to explicitly calculate the stochastic solution of Eq. (1) (see Appendix A)

$$\begin{aligned} Y(x)= \left\{ \begin{array}{ll} Y^{\textrm{I}}(x;E,W_0,W_1), &{}\quad 0<x \le l/2,\\ Y^{\textrm{I}}(x;E,W_0,W_1) +\dfrac{W_1-W_0}{384Ei}(l-2x)^4, &{}\quad l/2<x \le l, \end{array} \right. \end{aligned}$$

(6)

where

$$\begin{aligned} Y^{\textrm{I}}(x;E,W_0,W_1)=\dfrac{W_0 + 3 W_1}{16Ei}l^2 x^2-\dfrac{W_0 + W_1}{12Ei}l x^3+\dfrac{W_0}{24Ei}x^4. \end{aligned}$$

(7)

Observe that the solution Y(x) is differentiable, and therefore continuous at $x=l/2,$ which is interesting from an engineering point of view, since it informs us that the deflection varies smoothly at the spatial point, $x=l/2,$ where the load changes.

Since Y(x) is a piecewise function, the calculation of the 1-p.d.f. of Y(x) will be done separately in the two subdomains $0<x \le l/2$ and $l/2<x \le l$.

In the first case, we fix $0<x\le l/2$ and then we apply the RVT method, i.e. Theorem 1 to ${\textbf{U}}=\left( E, W_0, W_1 \right)$ to obtain the p.d.f. of the random vector ${\textbf{V}}=\left( V_1, V_2, V_3 \right) ,$ defined by the transformation ${\textbf{r}}:{\mathbb {R}}^3\rightarrow {\mathbb {R}}^3,$ whose components are conveniently defined as

$$\begin{aligned} v_1= & {} r_1(e,w_0,w_1) = Y^{\textrm{I}}(x;e,w_0,w_1), \\ v_2= & {} r_2(e,w_0,w_1) = w_0,\\ v_3= & {} r_3(e,w_0,w_1) = w_1. \end{aligned}$$

The inverse transformation of ${\textbf{r}},$ denoted by ${\textbf{s}}:{\mathbb {R}}^3\rightarrow {\mathbb {R}}^3,$ is given by

$$\begin{aligned} e= & {} s_1(v_1,v_2,v_3) = Z^{\textrm{I}}(x;v_1,v_2,v_3), \\ w_0= & {} s_2(v_1,v_2,v_3) = v_2,\\ w_1= & {} s_3(v_1,v_2,v_3) = v_3, \end{aligned}$$

where

$$\begin{aligned} Z^{\textrm{I}}(x;v_1,v_2,v_3)=\dfrac{1}{48v_1 i} \left( 3l^2 v_2 x^2 + 9 l^2 v_3 x^2 - 4 l v_2 x^3 - 4 l v_3 x^3 + 2v_2 x^4 \right) . \end{aligned}$$

(8)

To compute the p.d.f. of the random vector ${\textbf{V}},$ it is necessary to obtain the absolute value of the Jacobian of ${\textbf{s}},$

$$\begin{aligned} |{\mathcal {J}}| = \left| \frac{\partial s_1(v_1,v_2,v_3)}{\partial v_1}\right| = \left| \dfrac{1}{v_1}Z^{\textrm{I}}(x;v_1,v_2,v_3)\right| , \end{aligned}$$

which is well-defined and different from zero with probability 1 (w.p. 1) since $W_0,$ $W_1$ and E (hence $V_1,$ $V_2$ and $V_3$) are continuous random variables and $i \ne 0$. Applying the RVT method, stated in Theorem 1, we obtain the p.d.f. of the random vector ${\textbf{V}}=(V_1,V_2,V_3),$

$$\begin{aligned} f_{V_1,V_2,V_3} (v_1,v_2,v_3) = f_{E,W_0,W_1} (Z^{\textrm{I}}(x;v_1,v_2,v_3),v_2,v_3) \left| \dfrac{1}{v_1}Z^{\textrm{I}}(x;v_1,v_2,v_3)\right| . \end{aligned}$$

In the particular and important case that E,  $W_0$ and $W_1$ are independent random variables, the above expression writes

$$\begin{aligned} f_{V_1,V_2,V_3} (v_1,v_2,v_3) = f_{E}(Z^{\textrm{I}}(x;v_1,v_2,v_3)) f_{W_0}(v_2)f_{W_1} (v_3) \left| \dfrac{1}{v_1}Z^{\textrm{I}}(x;v_1,v_2,v_3)\right| . \end{aligned}$$

Notice that, according to the previous transformation ${\textbf{r}}:{\mathbb {R}}^3\rightarrow {\mathbb {R}}^3,$ the solution of model (1)–(2) corresponds to its first component $v_1,$ so marginalizing with respect to $V_2=W_0$ and $V_3=W_1,$ one obtains the 1-p.d.f. of Y(x), 

$$\begin{aligned} f_{Y(x)}(y)= & {} \int _{{\mathbb {R}}^2} f_{E}\left( \dfrac{1}{48y i} \left( 3l^2 w_0 x^2 + 9 l^2 w_1 x^2 - 4 l w_0 x^3 - 4 l w_1 x^3 + 2w_0 x^4 \right) \right) f_{W_0} (w_0)\nonumber \\{} & {} \cdot f_{W_1} (w_1) \dfrac{1}{48y^2 i} \left| \left( 3l^2 w_0 x^2 + 9 l^2 w_1 x^2 - 4 l w_0 x^3 - 4 l w_1 x^3 + 2w_0 x^4 \right) \right| {\textrm{d}} w_0\, {\textrm{d}} w_1. \end{aligned}$$

(9)

Remark 1

In the latter expression, it can be difficult to calculate the improper multi-integral using quadrature rules. For this reason and from a computational standpoint, it is interesting to express the above integral expression for the 1-p.d.f. in terms of the expectation operator w.r.t. random variables $(W_0,W_1)$ as follows

$$\begin{aligned} f_{Y(x)}(y)= & {} {\mathbb {E}}_{W_0,W_1}\left[ f_{E}\left( \dfrac{1}{48y i} \left( 3l^2 W_0 x^2 + 9 l^2 W_1 x^2 - 4 l W_0 x^3 - 4 l W_1 x^3 + 2W_0 x^4 \right) \right) \right. \nonumber \\{} & {} \left. \cdot \dfrac{1}{48y^2 i} \left| \left( 3l^2 W_0 x^2 + 9 l^2 W_1 x^2 - 4 l W_0 x^3 - 4 l W_1 x^3 + 2W_0 x^4 \right) \right| \right] , \quad 0<x \le l/2, \qquad \end{aligned}$$

(10)

since it permits applying Monte Carlo simulations to approximate the 1-p.d.f. We will use this remark throughout the article.

The computation of the 1-p.d.f. on the second piece of the domain $l/2<x \le l,$ can be calculated using a similar development as the one we have previously detailed step by step. We omit the technical details that lead to the following explicit expression in terms of the expectation operator

$$\begin{aligned} f_{Y(x)}(y)= & {} {\mathbb {E}}_{W_0,W_1}\left[ f_{E}\left( \dfrac{1}{384y i} \left( -l^4 W_0 +l^4 W_1 + 8 l^3 W_0 x -8 l^3 W_1 x + 96 l^2 W_1 x^2 \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. -64 l W_1 x^3 + 16 W_1 x^4 \right) \right) \cdot \dfrac{1}{384y^2 i} \left| \left( -l^4 W_0 +l^4 W_1 + 8 l^3 W_0 x -8 l^3 W_1 x \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + 96 l^2 W_1 x^2-64 l W_1 x^3 + 16 W_1 x^4 \right) \right| \right] , \quad l/2<x \le l. \end{aligned}$$

(11)

The maximum slope, S,  and the maximum deflection, D,  in a cantilever beam represent key information to account for safety and control measures (Öchsner 2021; Strømmen 2013; Mittelstedt 2021). In the case of a cantilever beam, it is clear they are calculated at the free end of the beam, i.e., $x=l$. It is easy to check that they are given by the following expressions

$$\begin{aligned} S = Y'(l) = \dfrac{l^3}{48E i}\left( W_0 + 7W_1 \right) , \quad D = Y(l) = \dfrac{l^4}{384Ei}\left( 7W_0 + 41W_1 \right) . \end{aligned}$$

(12)

In our context, S and D are random variables, whose respective p.d.f.’s, $f_{S}(s)$ and $f_{D}(d),$ provide important information about its probabilistic behaviour. As shown in Appendix B, the RVT technique permits explicitly calculating them

$$\begin{aligned} f_{S}(s) =\dfrac{l^3}{48 s^2 i} {\mathbb {E}}_{W_0,W_1}\left[ f_{E}\left( \dfrac{l^3}{48 s i} \left( W_0+7W_1 \right) \right) \left( W_0+7W_1 \right) \right] , \end{aligned}$$

(13)

and

$$\begin{aligned} f_{D}(d) = \dfrac{l^4}{384 d^2 i} {\mathbb {E}}_{W_0,W_1}\left[ f_{E}\left( \dfrac{l^4}{384 d i} \left( 7W_0+41W_1 \right) \right) \left( 7W_0+41W_1 \right) \right] . \end{aligned}$$

(14)

From an engineering point of view, it is useful to represent the response of the cantilever beam at each spatial point x in terms of static quantities as the bending moment, M(x),  and the shear force, V(x). These quantities can be obtained, respectively, in terms of the second and third order derivatives of the deflection (Strømmen 2020)

$$\begin{aligned} M(x) = -E i Y''(x), \end{aligned}$$

(15)

and

$$\begin{aligned} V(x) = -E i Y'''(x). \end{aligned}$$

(16)

In our setting problem, the following expressions for the bending moment and the shear force are, respectively, obtained

$$\begin{aligned} M(x)= & {} \left\{ \begin{array}{ll} -\frac{1}{8}(W_0 + 3 W_1) l^2 + \frac{1}{2}(W_0 + W_1) l x -\frac{1}{2} W_0 x^2, &{}\quad 0\le x \le l/2,\\ \left( -\frac{1}{2} l^2 + l x - \frac{1}{2} x^2\right) W_1, &{}\quad l/2<x \le l, \end{array} \right. \end{aligned}$$

(17)

$$\begin{aligned} V(x)= & {} \left\{ \begin{array}{ll} \frac{1}{2}(W_0 + W_1) l - W_0 x, &{}\quad 0\le x \le l/2,\\ \left( l - x \right) W_1, &{}\quad l/2<x \le l. \end{array} \right. \end{aligned}$$

(18)

We can observe that for $x=l$ both the bending moment and the shear force are null ($M(l)=0$ and $V(l)=0$). In addition, at $x=0$ they reach their maximum value (positive or negative).

Our goal again is to obtain the p.d.f. of (17) and (18). The RVT technique permits calculating the p.d.f. of the bending moment

$$\begin{aligned} f_{M(x)}(m)= \left\{ \begin{array}{ll} {\mathbb {E}}_{W_0}\left[ f_{W_1}\left( \dfrac{m-\left( -\frac{1}{8}l^2 + \frac{1}{2} l x - \frac{1}{2} x^2 \right) W_0}{\frac{1}{2} l x - \frac{3}{8} l^2} \right) \dfrac{1}{\left| \frac{1}{2} l x - \frac{3}{8} l^2 \right| } \right] , &{}\quad 0\le x \le l/2,\\ f_{W_1}\left( \dfrac{m}{ - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 } \right) \dfrac{1}{\left| - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 \right| }, &{}\quad l/2<x < l, \end{array} \right. \end{aligned}$$

(19)

and the p.d.f. of the shear force

$$\begin{aligned} f_{V(x)}(v)= \left\{ \begin{array}{ll} {\mathbb {E}}_{W_0}\left[ f_{W_1}\left( \dfrac{v-\left( \frac{1}{2} l - x \right) W_0}{\frac{1}{2} l} \right) \dfrac{2}{l} \right] , &{}\quad 0\le x \le l/2,\\ f_{W_1}\left( \dfrac{v}{l-x} \right) \dfrac{1}{l-x}, &{}\quad l/2<x < l. \end{array} \right. \end{aligned}$$

(20)

The development of expression (19) can be found in Appendix C. Expression (20) is obtained in a similar way.

Finally, we point out that the results derived throughout this section can be extended, using the same reasoning, to a finite number of loads including the case that the loads occupy different lengths on the beam.

3 Case II: Concentrated loads $P_j$ spanned on the beam

In this section we obtain the 1-p.d.f. of the solution stochastic process of model (1)–(2) in the important case that loads are applied at n different spatial points, $x_j, j=1,\ldots ,n,$ along the span of the beam. This case can be modelled taking

$$\begin{aligned} W(x) = \sum _{j=1}^{n}P_j \delta \left( x-x_j\right) , \end{aligned}$$

(21)

where $\delta (\cdot )$ denotes the Dirac delta distribution.

Let us assume that E and $P_j,$ $j=1,\dots ,n,$ are mutually independent random variables, being $f_E(e)$ and $f_{P_j}(p_j),$ $j=1,\ldots , n,$ their respective p.d.f.’s. In Fig. 2, we show a graphical scheme of the model. In a similar way that in the previous section, using the Laplace transform and its properties when computing the Laplace transform of the Dirac delta distribution, we can obtain the stochastic solution of model (1)–(2), via the following piecewise function

$$\begin{aligned} Y(x)= \left\{ \begin{array}{ll} Y^{\textrm{II}}(x;E,P_1,\dots ,P_n), &{} \quad 0< x\le x_1,\\ Y^{\textrm{II}}(x;E,P_1,\dots ,P_n) + \dfrac{1}{6Ei} \sum _{k=1}^{j}P_k \left( x-x_k \right) ^3, &{} \quad x_j<x\le x_{j+1},\\ &{}\quad j=1,\dots ,n-1,\\ Y^{\textrm{II}}(x;E,P_1,\dots ,P_n) + \dfrac{1}{6Ei} \sum _{k=1}^{n}P_k \left( x-x_k \right) ^3, &{} \quad x_n<x \le l, \end{array} \right. \end{aligned}$$

(22)

where

$$\begin{aligned} Y^{\textrm{II}}(x;E,P_1,\dots ,P_n)= \dfrac{1}{2Ei} x^2 \sum _{j=1}^{n}x_j P_j - \dfrac{1}{6Ei} x^3 \sum _{j=1}^{n}P_j. \end{aligned}$$

(23)

Fig. 2

Case II: Cantilever beam with concentrated random variables loads at different spatial points

Once we have obtained the stochastic solution, we will obtain the 1-p.d.f. of (22)–(23) taking advantage of the RVT method. As the solution is given by means of a piecewise function, the 1-p.d.f. will be determined by applying this method to each piece.

To this end, we first fix $0<x\le x_1$. Let us to define a mapping ${\textbf{r}}: {\mathbb {R}}^{n+1}\rightarrow {\mathbb {R}}^{n+1}$ that transforms the random vector ${\textbf{U}}=\left( E, P_1, \dots , P_n \right)$ into ${\textbf{V}}=\left( V_1, \dots , V_{n+1} \right)$ being

$$\begin{aligned} \begin{array}{ccccl} v_1 &{}=&{} r_1(e,p_1,\dots ,p_n) &{}=&{} Y^{\textrm{II}}(x;e,p_1,\dots ,p_n) , \\ v_2 &{}=&{} r_2(e,p_1,\dots ,p_n) &{}=&{} p_1,\\ \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots \\ v_{n+1} &{}=&{} r_{n+1}(e,p_1,\dots ,p_n) &{}=&{} p_{n}. \end{array} \end{aligned}$$

The inverse mapping of ${\textbf{r}}$ is given by the mapping ${\textbf{s}}: {\mathbb {R}}^{n+1}\rightarrow {\mathbb {R}}^{n+1}$ whose components are

$$\begin{aligned} \begin{array}{ccccl} e &{}=&{} s_1(v_1,v_2,\dots ,v_{n+1}) &{}=&{} Z^{\textrm{II}}(x;v_1,v_2,\dots ,v_{n+1}), \\ p_1 &{}=&{} s_2(v_1,v_2,\dots ,v_{n+1}) &{}=&{} v_2,\\ \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots \\ p_n &{}=&{} s_{n+1}(v_1,v_2,\dots ,v_{n+1}) &{}=&{} v_{n+1}, \end{array} \end{aligned}$$

where

$$\begin{aligned} Z^{\textrm{II}}(x;v_1,v_2,\dots ,v_{n+1})= -\dfrac{1}{6 v_1 i}\left( -3x^2 \sum _{j=1}^{n}x_j v_{j+1} +x^3\sum _{j=1}^{n}v_{j+1} \right) . \end{aligned}$$

Note that the absolute value of the Jacobian of ${\textbf{s}}$ is

$$\begin{aligned} |{\mathcal {J}}| = \left| \frac{\partial s_1(v_1,v_2,\dots ,v_{n+1})}{\partial v_1}\right| =\left| \frac{1}{v_1}Z^{\textrm{II}}(x;v_1,v_2,\dots ,v_{n+1})\right| \ne 0,\; {\mathrm {w.p.}}\, 1. \end{aligned}$$

Then, according to Theorem 1, the p.d.f. of ${\textbf{V}}$ is given by

$$\begin{aligned}{} & {} f_{V_1,\dots ,V_{n+1}} (v_1,\dots ,v_{n+1}) \\{} & {} \quad = f_{E}\left( Z^{\textrm{II}}(x;v_1,v_2,\dots ,v_{n+1})\right) f_{P_1}\left( v_2\right) \cdots f_{P_n}\left( v_{n+1}\right) \\{} & {} \qquad \cdot \left| \dfrac{1}{v_1}Z^{\textrm{II}}(x;v_1,v_2,\dots ,v_{n+1})\right| . \end{aligned}$$

Finally, the 1-p.d.f. of the stochastic solution (22)–(23), which is given by $V_1,$ is obtained by marginalizing w.r.t. $V_2 = P_1,$ $\dots ,$ $V_{n+1} = P_n,$

$$\begin{aligned} f_{Y(x)} (y)= & {} \int _{{\mathbb {R}}^n} f_{E}\left( -\dfrac{1}{6 y i} \left( -3 x^2 \sum _{j=1}^{n}x_j p_{j} + x^3\sum _{j=1}^{n}p_{j} \right) \right) f_{P_1}\left( p_1\right) \cdots f_{P_n}\left( p_{n}\right) \nonumber \\{} & {} \cdot \dfrac{1}{6 y^2 i} \left| -3x^2 \sum _{j=1}^{n}x_j p_j +x^3\sum _{j=1}^{n}p_j \right| {\textrm{d}} p_1 \cdots {\textrm{d}} p_n, \quad 0 < x\le x_1. \end{aligned}$$

(24)

The above expression involves a multidimensional integral that may be computationally unaffordable in practice. To overcome this drawback, it is interesting to notice that this expression can be rewritten using the operator expectation

$$\begin{aligned} f_{Y(x)} (y)= & {} {\mathbb {E}}_{P_1,\dots ,P_n}\left[ f_{E}\left( -\dfrac{1}{6 y i} \left( -3 x^2 \sum _{j=1}^{n}x_j P_{j} + x^3\sum _{j=1}^{n}P_{j} \right) \right) \right. \nonumber \\{} & {} \cdot \left. \dfrac{1}{6 y^2 i} \left| -3x^2 \sum _{j=1}^{n}x_j P_j +x^3\sum _{j=1}^{n}P_j \right| \right] , \quad 0 < x\le x_1. \end{aligned}$$

(25)

This expression is particularly useful when applying Monte Carlo simulations to calculate the 1-p.d.f.

The computation of the 1-p.d.f. on the other subdomains of the stochastic solution (see (22)) can be calculated similarly. For simplicity in the presentation, we here skip the technical details. Calculations lead to

$$\begin{aligned} f_{Y(x)} (y)= & {} {\mathbb {E}}_{P_1,\dots ,P_n}\left[ f_{E}\left( \dfrac{1}{6 y i} \left( 3 x^2 \sum _{k=1}^{n}x_k P_{k} - x^3\sum _{k=1}^{n}P_{k} + \sum _{k=1}^{j}\left( x-x_k \right) ^3 p_k\right) \right) \right. \nonumber \\{} & {} \cdot \left. \dfrac{1}{6 y^2 i} \left| 3x^2\! \sum _{k=1}^{n}x_k P_k -x^3\!\sum _{k=1}^{n}P_k +\sum _{k=1}^{j}\left( x-x_k \right) ^3 p_k \right| \right] , \nonumber \\{} & {} \qquad x_j<x\le x_{j+1},\; j=1,\dots ,n-1, \end{aligned}$$

(26)

and

$$\begin{aligned} f_{Y(x)} (y)= & {} {\mathbb {E}}_{P_1,\dots ,P_n}\left[ f_{E}\left( \dfrac{1}{6 y i} \left( 3 x^2 \sum _{k=1}^{n}x_k P_{k} - x^3\sum _{k=1}^{n}P_{k} + \sum _{k=1}^{n}\left( x-x_k \right) ^3 P_j\right) \right) \right. \nonumber \\{} & {} \cdot \left. \dfrac{1}{6 y^2 i} \left| 3x^2 \sum _{k=1}^{n}x_k P_k -x^3\sum _{k=1}^{n}P_k +\sum _{k=1}^{n}\left( x-x_k \right) ^3 P_k \right| \right] , \quad x_n<x \le l. \end{aligned}$$

(27)

Now, we are going to compute the p.d.f. of the slope, $S = Y'(l),$ and the maximum deflection, $D =Y(l),$ at free end as we did in the previous section. From (22)–(23) , it is easy to see that

$$\begin{aligned} S = \dfrac{1}{2Ei}\sum _{j=1}^{n}x_j^2 P_j, \quad D= \dfrac{1}{6Ei}\left( 3l\sum _{j=1}^{n}x_j^2 P_j - \sum _{j=1}^{n}x_j^3 P_j \right) , \end{aligned}$$

(28)

and applying the RVT technique with appropriate mappings, we can obtain the p.d.f.’s of S and D in terms of the expectation operator

$$\begin{aligned} f_{S}(s) = {\mathbb {E}}_{P_1,\dots ,P_n}\left[ f_{E}\left( \dfrac{1}{2 s i}\left( \sum _{j=1}^{n}x_j^2 P_j \right) \right) \dfrac{1}{2 s^2 i} \left| -\sum _{j=1}^{n}x_j^2 P_j \right| \right] , \end{aligned}$$

(29)

and

$$\begin{aligned} f_{D}(d) = {\mathbb {E}}_{P_1,\dots ,P_n}\left[ f_{E}\left( \dfrac{1}{6 d i}\left( 3 l \sum _{j=1}^{n}x_j^2 P_j -\sum _{j=1}^{n}x_j^3 P_j \right) \right) \dfrac{1}{6 d^2 i} \left| 3 l \sum _{j=1}^{n}x_j^2 P_j -\sum _{j=1}^{n}x_j^3 P_j \right| \right] . \end{aligned}$$

(30)

For the sake of completeness, we finish this section by calculating the p.d.f. of the bending moment and the shear force of the cantilever beam. Applying Eqs. (15) and (16) to (22)–(23), we obtain the expressions of the bending moment

$$\begin{aligned} M(x)= \left\{ \begin{array}{ll} \sum _{k=1}^{n}P_k (x - x_k), &{}\quad 0\le x \le x_1,\\ \sum _{k=j+1}^{n}P_k (x - x_k), &{}\quad x_j< x \le x_{j+1},\\ &{}\quad j=1,\dots , n-1,\\ 0, &{}\quad x_n < x \le l. \end{array} \right. \end{aligned}$$

(31)

and the shear force

$$\begin{aligned} V(x)= \left\{ \begin{array}{ll} \sum _{k=1}^{n}P_k, &{}\quad 0\le x \le x_1,\\ \sum _{k=j+1}^{n}P_k, &{}\quad x_j< x \le x_{j+1},\\ &{}\quad j=1,\dots , n-1,\\ 0, &{}\quad x_n < x \le l. \end{array} \right. \end{aligned}$$

(32)

Then, applying RVT technique to expressions (31) and (32) in a similar way to the previous section, we can obtain the p.d.f. of the bending moment

$$\begin{aligned} f_{M(x)}(m)= \left\{ \begin{array}{ll} {\mathbb {E}}_{P_1,\dots ,P_{n-1}}\!\left[ f_{P_n}\left( \frac{m-\sum _{k=1}^{n-1}P_k (x-x_k)}{x-x_n} \right) \dfrac{1}{\left| x-x_n \right| } \right] , &{}\quad 0\le x \le x_1,\\ {\mathbb {E}}_{ P_{j+1},\dots ,P_{n-1}}\left[ f_{P_n}\left( \frac{m-\sum _{k=j+1}^{n-1}P_k (x-x_k)}{x-x_n} \right) \dfrac{1}{\left| x-x_n \right| } \right] , &{}\quad x_j\le x \le x_{j+1},\\ &{}\quad j=1,\dots , n-2,\\ f_{P_n}\left( \dfrac{m}{x-x_n} \right) \dfrac{1}{\left| x-x_n \right| }, &{}\quad x_{n-1}< x \le x_n,\\ 0, &{}\quad x_n < x \le l, \end{array} \right. \end{aligned}$$

(33)

and the p.d.f. of the shear force

$$\begin{aligned} f_{V(x)}(v)= \left\{ \begin{array}{ll} {\mathbb {E}}_{P_2,\dots ,P_n}\left[ f_{P_1}\left( v - \sum _{k=2}^{n}P_k \right) \right] , &{}\quad 0\le x \le x_1,\\ {\mathbb {E}}_{P_{j+1},\dots ,P_{n-1}}\left[ f_{P_n}\left( v-\sum _{k=j+1}^{n-1}P_k \right) \right] , &{}\quad x_j< x \le x_{j+1}, \\ &{}\quad j=1,\dots , n-2,\\ f_{P_n}\left( v \right) , &{}\quad x_{n-1}< x \le x_n,\\ 0, &{}\quad x_n < x \le l. \end{array} \right. \end{aligned}$$

(34)

4 Case III: A load modelled via Brownian motion

In this section we assume that the density downward force, W(x),  is given by the sum of a deterministic value, $w_0$ (nominal value) and a certain random quantity that depends on the spatial point x on the beam. This latter quantity models the uncertainties due to heterogeneity of the beam. As for cantilever beams the key point is clearly located at the free end (in our case on the right end), then to mathematically analyse the problem we have chosen a stochastic process whose variance increases with its independent parameter (in our case the space x) and whose distribution is Gaussian. Specifically, we have selected a standard Wiener process (also termed Brownian motion), B(x),  to perform our mathematical analysis, although our subsequent approach also permits considering other stochastic processes to play the role of B(x). Therefore, we shall assume that

$$\begin{aligned} W(x)= w_0 + B(x), \quad 0<x\le l, \end{aligned}$$

(35)

being $B(x) = B(x, \omega )$ the standard Wiener process. Recall that, $\mu _B (x)={\mathbb {E}}[B(x)]=0$ and ${\mathbb {V}}[B(x)]=x,$ $\forall x >0,$ Kloeden and Platen (1992). Figure 3 shows a graphical representation of the problem.

Fig. 3

Case III: Cantilever beam with $w_0 + B(x)$ as varying load, being B(x) the Brownian motion

To perform the stochastic analysis by computing the 1-p.d.f. of the solution stochastic process in this case, we will take advantage of the Karhunen–Loève expansion of the Brownian motion (Lord et al. 2014, Ch. 5)

$$\begin{aligned} B(x) = \mu _B (x) + \sum _{j=1}^{\infty }\sqrt{\nu _j}\phi _j (x) \xi _j (\omega ),\;\; \omega \in \Omega , \, \, 0<x \le l, \end{aligned}$$

(36)

where $\mu _B (x)=0,$ $\xi _j (\omega )$ are independent and identically distributed standard Gaussian random variables, $\xi _j(\omega )\sim {\textrm{N}}(0,1),$ $j=1, 2, \dots ,$ and $\{\nu _j,\phi _j(x)\}$ are the eigenpairs (eigenvalues and eigenfunctions) obtained when solving the following homogeneous Fredholm integral equation of second kind

$$\begin{aligned} \int _{0}^{l} c_B(x_1,x_2)\, \phi _j(x_2)\, {\textrm{d}} x_2 = \nu _j\,\phi _j(x_1). \end{aligned}$$

Here, $c_B(x_1,x_2)$ is the covariance function of the Brownian motion, B(x). It can be seen that (Lord et al. 2014, p. 216)

$$ c_B(x_1,x_2)=\min (x_1,x_2), \qquad (x_1,x_2)\in [0,l]\times [0,l], $$

and

$$\begin{aligned} \nu _j = \dfrac{4 l^2}{\pi ^2(2j-1)^2}, \quad \phi _j(x)=\sqrt{\frac{2}{l}}\sin {\left( \frac{(2j-1)\pi }{2l}x \right) }, \quad j=1, 2, \dots \end{aligned}$$

To obtain an approximation of the 1-p.d.f. of the solution stochastic process of model (1)–(2) with W(x) given by (35), we will consider the approximation of B(x) obtained by truncating its Karhunen–Loève expansion at a finite order, say N. Then, the model is approximated via the differential equation

$$\begin{aligned} \frac{{\textrm{d}}^4 Y(x)}{{\textrm{d}} x^4}=\frac{1}{Ei}\left( w_0 + \sum _{j=1}^{N}\sqrt{\nu _j}\phi _j (x) \xi _j (\omega ) \right) , \quad 0<x<l. \end{aligned}$$

(37)

The solution stochastic process of model (37) together with the boundary conditions (2) can be obtained using the Laplace transform. After somewhat technical computations, we obtain

$$\begin{aligned} Y(x)= & {} \dfrac{1}{Ei}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}-\sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) \right. \nonumber \\{} & {} \left. + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) + \dfrac{w_0}{24}x^4 \right. \nonumber \\{} & {} \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) , \quad 0<x \le l, \end{aligned}$$

(38)

where $b_j = \dfrac{(2j-1)\pi }{2l}$.

We fix $0<x<l$ and we assume that E and $\xi _j,$ $j=1\dots ,N$ are independent continuous random variables with p.d.f.’s given by $f_E(e),$ $f_{\xi _j}(\xi _j),$ $j=1\dots ,N,$ respectively. Now, we will apply the RVT method taking in Theorem 1 as ${\textbf{U}}=\left( E, \xi _1,\dots , \xi _N \right)$ to obtain the p.d.f. of the random vector ${\textbf{V}}=\left( V_1, V_2, \dots , V_{N+1} \right) ,$ defined by the transformation ${\textbf{r}}:{\mathbb {R}}^{N+1}\rightarrow {\mathbb {R}}^{N+1},$ whose components are defined by

$$\begin{aligned} v_1= & {} r_1(e,\xi _1,\dots ,\xi _N) = \dfrac{1}{e i}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right. \right. \nonumber \\{} & {} \left. \left. -\sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) \right. \nonumber \\{} & {} \left. + \dfrac{w_0}{24}x^4 + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) , \nonumber \\ v_2= & {} r_2(e,\xi _1,\dots ,\xi _N) = \xi _1,\nonumber \\ \vdots \nonumber \\ v_{N+1}= & {} r_{N+1}(e,\xi _1,\dots ,\xi _N) = \xi _{N}. \end{aligned}$$

The inverse transformation of ${\textbf{r}},$ ${\textbf{s}}:{\mathbb {R}}^{N+1}\rightarrow {\mathbb {R}}^{N+1},$ is given by

$$\begin{aligned} e= & {} s_1(v_1,v_2,\dots ,v_{N+1}) = \dfrac{1}{v_1 i}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}v_{j+1} \dfrac{1-\cos {(b_j l)}}{b_j^2} \right. \right. \nonumber \\{} & {} \left. \left. -\sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) \right. \nonumber \\{} & {} \left. + \dfrac{w_0}{24}x^4 + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) , \nonumber \\ \xi _1= & {} s_2(v_1,v_2,\dots ,v_{N+1}) = v_2,\\ \vdots \nonumber \\ \xi _N= & {} s_{N+1}(v_1,v_2,\dots ,v_{N+1}) = v_{N+1}. \end{aligned}$$

The absolute value of the Jacobian of the inverse mapping ${\textbf{s}}$ writes

$$\begin{aligned} |{\mathcal {J}}|= & {} \left| \frac{\partial s_1(v_1,v_2,\dots ,v_{N+1})}{\partial v_1}\right| = \left| -\dfrac{1}{v_1^2 i}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}v_{j+1} \dfrac{1-\cos {(b_j l)}}{b_j^2}\right. \right. \right. \\{} & {} \left. \left. \left. -\sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) \right. \right. \\{} & {} \left. \left. + \dfrac{w_0}{24}x^4 + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}v_{j+1} \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) \right| \ne 0,\; {\mathrm {w.p.}}\, 1. \end{aligned}$$

Similarly as developed in cases I and II, the following expression for the 1-p.d.f. of the solution stochastic process is obtained

$$\begin{aligned} f_{Y(x)}(y)= & {} {\mathbb {E}}_{\xi _1,\dots , \xi _{N}}\left[ f_{E}\left( \dfrac{1}{E i}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right. \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. \left. -\sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + \dfrac{w_0}{24}x^4 + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) \right) \right. \nonumber \\{} & {} \left. \cdot \left| -\dfrac{1}{E^2 i}\left( \dfrac{x^2}{2} \left( \dfrac{w_0}{2}l^2+l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}-\sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3}\right) \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + \dfrac{x^3}{6} \left( -w_0 l - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} \right) + \dfrac{w_0}{24}x^4 \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j x+b_j^3 x^3+6\sin {(b_j x)}}{6 b_j^5} \right) \right| \right] , \quad 0<x\le l. \end{aligned}$$

(39)

To complete this case III with same information as the one presented in cases I and II, now, we compute the p.d.f.’s of the maximum slope and deflection at the free end of the beam and the p.d.f.’s of the bending moment and the shear force. For the sake of simplicity, we directly present the obtained results. For the maximum slope

$$\begin{aligned} S =Y'(l)= & {} \dfrac{1}{Ei}\left( \dfrac{w_0}{8}l^3 + \dfrac{1}{2}l^2 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}-l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} \right. \nonumber \\{} & {} \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-2+b_j^2 l^2+2\cos {(b_j l)}}{2 b_j^4} \right) , \end{aligned}$$

(40)

its p.d.f. is given by

$$\begin{aligned} f_{S}(s)= & {} {\mathbb {E}}_{\xi _1,\dots ,\xi _N}\left[ f_{E}\left( \dfrac{1}{s i}\left( \dfrac{w_0}{8}l^3 + \dfrac{1}{2}l^2 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}\right. \right. \right. \nonumber \\{} & {} -l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} \nonumber \\{} & {} \left. \left. \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-2+b_j^2 l^2+2\cos {(b_j l)}}{2 b_j^4} \right) \right) \right. \left| -\dfrac{1}{\theta ^2 i}\left( \dfrac{w_0}{8}l^3 + \dfrac{1}{2}l^2 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}\right. \right. \nonumber \\{} & {} \left. \left. \left. -l\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-2+b_j^2 l^2+2\cos {(b_j l)}}{2 b_j^4} \right) \right| \right] , \end{aligned}$$

(41)

and, for the maximum deflection

$$\begin{aligned} D= Y(l)= & {} \dfrac{1}{Ei}\left( \dfrac{w_0}{8}l^4 + \dfrac{1}{3}l^3 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}-\dfrac{1}{2}l^2\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} \right. \nonumber \\{} & {} \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j l+b_j^3 l^3+6\sin {(b_j l)}}{6 b_j^5} \right) , \end{aligned}$$

(42)

its p.d.f. is given by

$$\begin{aligned} f_{D}(d)= & {} {\mathbb {E}}_{\xi _1,\dots ,\xi _N}\left[ f_{E}\left( \dfrac{1}{d i}\left( \dfrac{w_0}{8}l^4 + \dfrac{1}{3}l^3 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}\right. \right. \right. \nonumber \\{} & {} -\frac{1}{2}l^2\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} \nonumber \\{} & {} \left. \left. + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j l+b_j^3 l^3+6\sin {(b_j l)}}{6 b_j^5} \right) \right) \left| -\dfrac{1}{\delta ^2 i}\left( \dfrac{w_0}{8}l^4 + \dfrac{1}{3}l^3 \sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2}\right. \right. \nonumber \\{} & {} \left. \left. \left. -\frac{1}{2}l^2\sqrt{\frac{2}{l}}\sum _{j=1}^{N}\xi _j \dfrac{b_j l-\sin {(b_j l)}}{b_j^3} + \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{-6b_j l+b_j^3 l^3+6\sin {(b_j l)}}{6 b_j^5} \right) \right| \right] . \end{aligned}$$

(43)

Applying Eqs. (15) and (16) to (38), the following expressions are, respectively, obtained for the bending moment

$$\begin{aligned} M(x)= & {} -\frac{w_0}{2}\left( l^2 + x^2 - 2 x l \right) - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{1-\cos {(b_j l)}}{b_j^2} (l-x)\nonumber \\{} & {} - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{b_j (x-l) + \sin {(b_j l)} - \sin {(b_j x)}}{b_j^3}, \end{aligned}$$

(44)

and the shear force

$$\begin{aligned} V(x) = -w_0 ( x - l) - \sqrt{\frac{2}{l}} \sum _{j=1}^{N}\xi _j \dfrac{\cos {(b_j l)}-\cos {(b_j x)}}{b_j^2}. \end{aligned}$$

(45)

Its p.d.f.’s are given by

$$\begin{aligned} f_{M(x)}(m)= & {} {\mathbb {E}}_{\xi _2,\dots ,\xi _{N}}\left[ f_{\xi _1}\left( \left( m + \frac{w_0}{2}\left( l^2 + x^2 - 2 x l \right) \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + \sqrt{\frac{2}{l}} \sum _{j=2}^{N}\xi _j \left( \dfrac{1-\cos (b_j l)}{b_j^2}(l-x) + \dfrac{b_j (x-l) + \sin (b_j l)-\sin (b_j x)}{b_j^3} \right) \right) \right. \right. \nonumber \\{} & {} \left. \left. \cdot \left( -\sqrt{\frac{2}{l}} \left( \dfrac{1-\cos (b_1 l)}{b_1^2}(l-x) + \dfrac{b_1 (x-l) + \sin (b_1 l)-\sin (b_1 x)}{b_1^3}\right) \right) ^{-1} \right) \right. \nonumber \\{} & {} \left. \cdot \left| -\sqrt{\frac{2}{l}} \left( \dfrac{1-\cos (b_1 l)}{b_1^2}(l-x) + \dfrac{b_1 (x-l) + \sin (b_1 l)-\sin (b_1 x)}{b_1^3}\right) \right| ^{-1} \right] ,\nonumber \\ \end{aligned}$$

(46)

and

$$\begin{aligned} f_{V(x)}(v)= & {} {\mathbb {E}}_{\xi _2,\dots ,\xi _{N}}\left[ f_{\xi _1}\left( \left( v + w_0(x - l) + \sqrt{\frac{2}{l}} \sum _{j=2}^{N}\xi _j \dfrac{\cos (b_j l)-\cos (b_j x)}{b_j^2} \right) \right. \right. \nonumber \\{} & {} \left. \left. \cdot \left( -\sqrt{\frac{2}{l}} \left( \dfrac{\cos (b_1 l)-\cos (b_1 x)}{b_1^2} \right) \right) ^{-1} \right) \left| -\sqrt{\frac{2}{l}} \left( \dfrac{\cos (b_1 l)-\cos (b_1 x)}{b_1^2} \right) \right| ^{-1} \right] , \end{aligned}$$

(47)

respectively.

Remark 2

In the previous analysis we have calculated the 1-p.d.f. of several quantities of interest (deflection, shear force and bending moment) of the beam thanks to the application of the RVT method. It is worthwhile to point out that this technique could also be applied to determine the second probability density function (2-p.d.f.) of the above mentioned quantities (Soong 1973). This would allow us to quantify further statistical properties, as for instance, the correlation of the deflection at two different spatial points. Here, we omit this analysis since is very similar but involving cumbersome mathematical expressions.

5 Numerical examples

This section is devoted to illustrate the above theoretical findings. We show three different examples corresponding to the theoretical results obtained in cases I–III developed in each one of the previous Sects. 2–4, respectively. We take the following data for the deterministic parameters of the model (1): $l = 10~\text {m}$ and $i = 722~\text {cm}^4,$ while for the random parameters, we will assume that the Young’s modulus of elasticity, E,  has a truncated Gaussian distribution with mean $\mu _E = 210\times 10^9$ ${\textrm{Pa}}$ (Pascals) and variance $\sigma ^2_E=420\times 10^7$ $\text {Pa}^2,$ i.e. $E \sim \text {N}_{{\mathcal {T}}}(\mu _E; \sigma ^2_E)$ where ${\mathcal {T}}=[\mu _E-k \sigma _E,\mu _E+k \sigma _E] = [209.9993\times 10^9, 210.0006\times 10^9],$ taking $k=10$. The particular form of stochastic process W(x) will be specified below in each example.

Example 1

(Case I) This example corresponds to Sect. 2, where the cantilever beam supports two different loads. We assume that the loads are defined by random variables whose distributions are Gaussian, $W_0\sim {\textrm{N}}\left( \mu _{W_0}=40; \sigma _{W_0}^2=0.4 \right)$ and $W_1\sim {\textrm{N}}\left( \mu _{W_1}=20; \sigma _{W_1}^2=0.2 \right)$.

Figure 4 shows the graphical representation of the 1-p.d.f. given by expressions (10) and (11) at different spatial points $x \in \{1,\ldots ,10\}$. As it is expected, we can observe that the variance increases as x does.

Fig. 4

1-p.d.f., $f_{Y(x)}(y),$ of the solution stochastic process (6), computed by (10) and (11), at different spatial points $x \in \{1,\ldots ,10\}$. Example 1

In Fig. 5, we show the plots of the p.d.f.’s of the maximum slope, $f_{S}(s)$ and the maximum deflection, $f_{D}(d),$ at free end of the cantilever beam. Computations have been carried out by expressions (13) and (14), respectively. In Table 1, we include the numerical results for the mean and the standard deviation of random variables S and D.

Fig. 5

Left: P.d.f. of the maximum slope, $f_{S}(s),$ at free end using expression (13). Right: P.d.f. of the maximum deflection, $f_{D}(d),$ at free end using (14). Example 1

Table 1 Mean and standard deviation of the maximum slope and the maximum deflection of Y(x) at the free end of the beam

In Figs. 6 and 7 we show a graphical representation of the p.d.f.’s of the bending moment, $f_{M(x)}(m)$ and the shear force, $f_{V(x)}(v),$ respectively, at different spatial points $x \in \{0,\ldots ,9\}$. Recall that at the end of the beam ($l=10$), $M(l)=0$ and $V(l)=0$. We can observe in both figures that the variance decreases as the position increases.

Fig. 6

P.d.f. of the bending moment, $f_{M(x)}(m),$ using expression (19) at different spatial points $x \in \{0,\ldots ,9\}$. Example 1

Fig. 7

P.d.f. of the shear force, $f_{V(x)}(v),$ using expression (20) at different spatial points $x \in \{0,\ldots ,9\}$. Example 1

Finally, in Fig. 8 we show a graphical representation of the mean and standard deviation functions of the solution stochastic process, Y(x). They have been calculated by expressions (3) and (4), where $f_{Y(x)}(y)$ is given by (10)–(11).

Fig. 8

Mean ($\mu$) plus/minus 2 standard deviations ($\sigma$) of the solution stochastic process, Y(x). Example 1

Example 2

(Case II) Here we illustrate the theoretical results for a cantilever beam subject to concentrated loads at specific spatial points. More specifically, we assume that four random loads are located at $x_j = 2,5,7,9,$ $j=1,2,3,4,$ respectively. We will assume that all the loads are characterized by a common Gaussian distribution, $P_j \sim {\textrm{N}}(\mu _{P_j}=20; \sigma _{P_j}^2= 0.02),$ $j=1,\ldots ,4$.

In Fig. 9 we have graphically represented the 1-p.d.f., $f_{Y(x)}(y),$ computed by (25)–(27), at the spatial positions $x=j,$ $j=1,\ldots ,10$.

Fig. 9

1-p.d.f., $f_{Y(x)}(y),$ of the solution stochastic process (22), computed by (25)–(27), at the different spatial points $x= j,$ $j=1,2,\ldots ,10$. Example 2

In Fig. 10, we have plotted the p.d.f. of the maximum slope, $f_{S}(s),$ and of the maximum deflection, $f_{D}(d),$ at the end of the cantilever beam. These two densities are given by (29) and (30), respectively. In Table 2, we present the numerical results of the mean and standard deviation of these two random variables.

Fig. 10

Left: P.d.f., $f_{S}(s),$ of the maximum slope at free end. Right: P.d.f., $f_{D}(d),$ of the maximum deflection at the free end. Example 2

Table 2 Mean and standard deviation of the maximum slope and the maximum deflection of Y(x)

In Fig. 11 we show a graphical representation of the p.d.f. of the bending moment, $f_{M(x)}(m)$ at different spatial points $x \in \{0,\ldots ,8\}$. Due to expression (31) at instants $x \in \{9, 10\},$ the value of the p.d.f. computed by (33) is zero. Figure 12 shows the plot of the p.d.f.’s of the shear force, $f_{V(x)}(v)$. Notice in this graphical representation that several p.d.f.’s match, namely, $f_{V(0)}(v)=f_{V(1)}(v)=f_{V(2)}(v),$ $f_{V(3)}(v)=f_{V(4)}(v)=f_{V(5)}(v),$ $f_{V(6)}(v)=f_{V(7)}(v)$ and $f_{V(8)}(v)=f_{V(9)}(v)$. This is a consequence of the definition of the shear force (see Eq. (32)) and the spatial points, $x_j = 2,5,7,9$ (meters), $j=1,2,3,4,$ where the loads have been placed. We can observe in Figs. 11 and 12 that the variance decreases as the position increases.

Fig. 11

P.d.f. of the bending moment, $f_{M(x)}(m),$ using expression (33) at different spatial points $x \in \{0,\ldots ,8\}$. Example 2

Fig. 12

P.d.f. of the shear force, $f_{V(x)}(v),$ using expression (34) at different spatial points $x \in \{0,\ldots ,9\}$. Example 2

Finally, in Fig. 13 we show the mean plus/minus 2 standard deviations of the solution stochastic process, Y(x).

Fig. 13

Mean ($\mu$) plus/minus 2 standard deviations ($\sigma$) of the solution stochastic process, Y(x). Example 2

Example 3

(Case III) To illustrate the findings obtained in Sect. 4, we will consider that the density of downward force acting perpendicularly on the beam of length $l=10\,\textrm{m}$ is given by $W(x) = w_0+ B(x),$ $0<x<l,$ $w_0=20,$ and we will consider a Karhunen–Loève expansion truncated at order N to approximate the Brownian motion, B(x).

In Fig. 14, we show a graphical representation of the 1-p.d.f., $f_{Y(x)}(x),$ given by (39) for different spatial points $x \in \{1,\ldots ,10\}$ considering as truncating order $N=1$ (later we justify this approximation is enough to achieve reliable results). We can notice that the higher the position, the higher the variability as expected, since the variability of the Brownian motion, B(x),  increases with x.

Fig. 14

1-p.d.f., $f_{Y(x)}(y),$ of the solution stochastic process (38), computed by (39), at different spatial position $x \in \{1,\ldots ,10\}$ of the cantilever beam considering an approximation of the Brownian motion, B(x),  via a Karhunen–Loève expansion truncated at order $N=1$. Example 3

In Fig. 15, we show the p.d.f. of the maximum slope and maximum deflection at the end of the cantilever beam. From both p.d.f.’s, we have calculated, in Table 3, approximations of the mean and standard deviation of these two random variables.

Fig. 15

Left: P.d.f., $f_{S}(s),$ of the maximum slope at free end. Right: P.d.f., $f_{D}(d),$ of the maximum deflection at the free end. Example 3

In Figs. 16 and 17, we show, respectively, a graphical representation of the p.d.f.’s of the bending moment, $f_{M(x)}(m),$ and the p.d.f.’s of the shear force, $f_{V(x)}(v),$ at different spatial points $x \in \{0,\ldots ,9\}$. Recall that at $x=0$ reach their maximum value (positive or negative) and at $x=l=10,$ $M(l)=0$ and $V(l)=0$. To compute these p.d.f.’s we have considered again the truncating order $N=1$. This causes expressions (46) and (47) to change their structure a slight bit. So, the p.d.f. of the bending moment for $N=1$ is given by

$$\begin{aligned} f_{M(x)}(m)= & {} f_{\xi _1}\left( \left( m + \frac{w_0}{2}\left( l^2 + x^2 - 2 x l \right) \right) \left( -\sqrt{\frac{2}{l}} \left( \dfrac{1-\cos (b_1 l)}{b_1^2}(l-x) \right. \right. \right. \nonumber \\{} & {} \left. \left. \left. + \dfrac{b_1 (x-l) + \sin (b_1 l)-\sin (b_1 x)}{b_1^3}\right) \right) ^{-1} \right) \nonumber \\{} & {} \cdot \left| -\sqrt{\frac{2}{l}} \left( \dfrac{1-\cos (b_1 l)}{b_1^2}(l-x) + \dfrac{b_1 (x-l) + \sin (b_1 l)-\sin (b_1 x)}{b_1^3}\right) \right| ^{-1}, \nonumber \\ \end{aligned}$$

(48)

and the p.d.f. of the shear force is given by

$$\begin{aligned} f_{V(x)}(v)= & {} f_{\xi _1}\left( \left( v + w_0(x - l) \right) \left( -\sqrt{\frac{2}{l}} \left( \dfrac{\cos (b_1 l)-\cos (b_1 x)}{b_1^2} \right) \right) ^{-1} \right) \nonumber \\{} & {} \cdot \left| -\sqrt{\frac{2}{l}} \left( \dfrac{\cos (b_1 l)-\cos (b_1 x)}{b_1^2} \right) \right| ^{-1}. \end{aligned}$$

(49)

Table 3 Mean and standard deviation of the maximum slope and deflection at free end of the cantilever beam, obtained via the p.d.f. (41) and (43), respectively

Again, we can observe in both figures that the variance decreases as the position increases.

Fig. 16

P.d.f. of the bending moment, $f_{M(x)}(m),$ using expression (46) at different spatial points $x \in \{0,\ldots ,8\}$. Example 3

Fig. 17

P.d.f. of the shear force, $f_{V(x)}(v),$ using expression (47) at different spatial points $x \in \{0,\ldots ,9\}$. Example 3

In Fig. 18, we show the mean ($\mu$) plus/minus 2 standard deviations of the solution stochastic process, Y(x),  on the whole spatial domain.

Fig. 18

Mean ($\mu$) plus/minus 2 standard deviations of the solution stochastic process, Y(x). Example 3

Finally, in Table 4, we show a comparison of the values of the mean and standard deviation of the maximum deflection of the cantilever beam, D,  considering different orders of truncation, $N \in \{1,2,3,10,50\},$ of B(x) via a Karhunen–Loève expansion. We can observe that approximations are very similar with $N=1$. This justifies that our previous calculations have been carried out computations with this order. Similar conclusions are derived from the p.d.f.’s as we can observe from Fig. 19.

Table 4 Mean and standard deviation of the maximum deflection, D,  at free end, obtained via the p.d.f. (43) for different values of the truncation order, N,  to approximate the Brownian motion by its Karhunen–Loève expansion, $N \in \{1,2,3,10,50\}$

Fig. 19

P.d.f. of the maximum deflection at free end, $f_{D}(d),$ for different values of the truncation order, N,  to approximate the Brownian motion by its Karhunen–Loève expansion $N\in \{1,2,3,10,50\}$. Example 3

6 Conclusions

Throughout the paper, we have studied, from a probabilistic standpoint, a foundational model to describe the deflection of a cantilever beam subject to different types of loads, which are relevant in the civil engineering literature. Our analysis has several advantages, first, it permits computing not only the mean and the standard deviation of deflection, but its probability density function that provides a fuller description. Secondly, our theoretical findings have been obtained under very general hypotheses since we have considered a complete randomization of model parameters (the density of downward force acting perpendicularly on the beam and the Young’s modulus, denoted by W(x) and E,  respectively). Even more, the results obtained in every of the three cases analysed in the paper have been established assuming arbitrary probability density functions for the random parameters involved in the model. This gives a great generality to our results. Besides, we have probabilistically characterized, via the corresponding densities, important quantities of interest as the slope, maximum deflection, bending moment and shear force of the cantilever beam under mild hypotheses. Furthermore, we have shown that the Random Variable Transformation technique provides a comprehensive, systematic and unifying tool to obtain the mathematical results in a variety of scenarios, allowing us to obtain general formulas that are very useful to carry out computations as shown in the examples. We think that this approach can be promising to deal with, in future works, the stochastic analysis of the deflection of other kinds of beams than cantilevers.

Appendices

Appendix A: Solution stochastic process: case I

As we have seen in Sect. 2, to obtain the 1-p.d.f., $f_{Y(x)}(y),$ we need to explicitly calculate the stochastic solution of Eq. (1)–(2). To this end several techniques can be applied, here we shall apply the Laplace transform because of the well-known advantages of this integral transform to compactly obtain the solution Y(x) when W(x) is defined via a piecewise function representing the loads spanned on the beam (Iwiński 1967). Furthermore, the resulting expression of the solution is particularly useful to apply the RVT technique, as it will be apparent throughout our subsequent development. To apply the Laplace transform, we first need to extend the definition of (5) as

$$\begin{aligned} W(x)= \left\{ \begin{array}{ll} W_0, &{}\quad 0<x \le l/2,\\ W_1, &{}\quad l/2<x \le l,\\ 0, &{}\quad x>l. \end{array} \right. \end{aligned}$$

The above expression can be written in terms of a unitary Heaviside function, ${\mathcal {U}}(x),$

$$\begin{aligned} W(x) = W_0 \left( {\mathcal {U}}(x)-{\mathcal {U}}(x-l/2) \right) + W_1 \left( {\mathcal {U}}(x-l/2) - {\mathcal {U}}(x-l) \right) . \end{aligned}$$

(A1)

Now, applying the Laplace transform, $y(s)={\mathcal {L}}\{Y(x)\}(s)=\int _{0}^{\infty } \textrm{e}^{-s x} Y(x) \, {\textrm{d}}x$ to (1)–(2) with W(x) given by (A1), using the well-known properties of this integral transform (Schiff 1999, Ch. 2) and isolating y(s),  one obtains

$$\begin{aligned} y(s) = \dfrac{1}{s^4}\left( s c_1 + c_2 +\dfrac{1}{E i} \left( \dfrac{W_0}{s} - \dfrac{W_0}{s} \text {e}^{-\frac{l s}{2}} + \dfrac{W_1}{s} \text {e}^{-\frac{l s}{2}} - \dfrac{W_1}{s} \text {e}^{-l s} \right) \right) , \end{aligned}$$

(A2)

where $c_1 = Y''(0)$ and $c_2 = Y'''(0)$. The next step is to apply the inverse Laplace transform to (A2)

$$\begin{aligned} {\mathcal {L}}^{-1}\{y(s)\} = {\mathcal {L}}^{-1}\left\{ \dfrac{1}{s^3}c_1 + \dfrac{1}{s^4}c_2 + \dfrac{1}{E i} \left( W_0\dfrac{1}{s^5} + \left( W_1 - W_0 \right) \dfrac{1}{s^5}\text {e}^{-\frac{l s}{2}} - W_1 \dfrac{1}{s^5}\text {e}^{-l s} \right) \right\} . \end{aligned}$$

Then, using the inverse Laplace properties (Schiff 1999), and applying the boundary conditions $Y''(l)=0$ and $Y'''(l)=0$ to obtain the values of $c_1$ and $c_2,$ we obtain the solution stochastic process of model (1)–(2)

$$\begin{aligned} Y(x)= \left\{ \begin{array}{lc} Y^{\textrm{I}}(x;E,W_0,W_1), &{}\quad 0<x \le l/2,\\ Y^{\textrm{I}}(x;E,W_0,W_1) +\dfrac{W_1-W_0}{384Ei}(l-2x)^4, &{}\quad l/2<x \le l, \end{array} \right. \end{aligned}$$

where

$$\begin{aligned} Y^{\textrm{I}}(x;E,W_0,W_1)=\dfrac{l^2 W_0 + 3l^2 W_1}{16Ei}x^2-\dfrac{l W_0 + l W_1}{12Ei}x^3+\dfrac{W_0}{24Ei}x^4. \end{aligned}$$

Appendix B: P.d.f. of the maximum slope at free end: case I

As we have seen in Sect. 2, the expression of the maximum slope at the free end in the cantilever beam is given by

$$\begin{aligned} S = \dfrac{l^3}{48E i}\left( W_0 + 7W_1 \right) , \end{aligned}$$

(B3)

where $W_0,$ $W_1$ and E are independent continuous random variables whose p.d.f.’s are denoted by $f_E(e),$ $f_{W_0}(w_0)$ and $f_{W_1}(w_1) ,$ respectively.

To obtain the p.d.f. of S,  we apply the RVT method stated in Theorem 1 taking ${\textbf{U}}=\left( E, W_0, W_1 \right)$ to obtain the p.d.f. of the random vector ${\textbf{V}}=\left( V_1, V_2, V_3 \right) ,$ defined by the transformation ${\textbf{r}}:{\mathbb {R}}^3\rightarrow {\mathbb {R}}^3,$ whose components are given by

$$\begin{aligned} v_1= & {} r_1(e,w_0,w_1) = \frac{l^3}{48 e i}\left( w_0 + 7w_1 \right) , \\ v_2= & {} r_2(e,w_0,w_1) = w_0,\\ v_3= & {} r_3(e,w_0,w_1) = w_1. \end{aligned}$$

The inverse transformation of ${\textbf{r}},$ ${\textbf{s}}: {\mathbb {R}}^{3}\rightarrow {\mathbb {R}}^{3},$ writes

$$\begin{aligned} e= & {} s_1(v_1,v_2,v_3) = \frac{l^3}{48 v_1 i}\left( v_2 + 7v_3 \right) , \\ w_0= & {} s_2(v_1,v_2,v_3) = v_2,\\ w_1= & {} s_3(v_1,v_2,v_3) = v_3. \end{aligned}$$

The absolute value of the Jacobian of the inverse transformation is

$$\begin{aligned} | {\mathcal {J}}| = \left| \frac{\partial s_1(v_1,v_2,v_3)}{\partial v_1}\right| = \dfrac{l^3}{48 v^2_1 i}\left( v_2 + 7v_3 \right) \ne 0,\quad {\mathrm {w.p.}} \, 1. \end{aligned}$$

So, the p.d.f. of ${\textbf{V}}$ is given by

$$\begin{aligned} f_{V_1,V_2,V_3} (v_1,v_2,v_3) = f_{E} \left( \frac{l^3}{48 v_1 i}\left( v_2 + 7v_3 \right) \right) f_{W_0}(v_2)f_{W_1}(v_3) \frac{l^3}{48 v^2_1 i}\left( v_2 + 7v_3 \right) . \end{aligned}$$

To obtain the p.d.f. of (B3), which is given by $V_1,$ we marginalize w.r.t. $V_2 = W_0$ and $V_{3} = W_1,$ turning out

$$\begin{aligned} f_{S}(s) =\dfrac{l^3}{48 s^2 i} \int _{{\mathbb {R}}^2} f_{E}\left( \dfrac{l^3}{48 s i} \left( w_0 + 7 w_1 \right) \right) f_{W_0} (w_0)f_{W_1} (w_1) \left( w_0 + 7 w_1 \right) {\textrm{d}} w_0\, {\textrm{d}} w_1. \end{aligned}$$

This p.d.f. can be expressed by the expectation of the random vector $(W_0,W_1)$ as

$$\begin{aligned} f_{S}(s) = \dfrac{l^3}{48 s^2 i} {\mathbb {E}}_{W_0,W_1}\left[ f_{E}\left( \dfrac{l^3}{48 s i} \left( W_0 + 7 W_1 \right) \right) \left( W_0+7W_1\right) \right] . \end{aligned}$$

Appendix C: P.d.f. of the bending moment: case I

As we have seen in Sect. 2, the bending moment in the cantilever beam subject to two random loads is given by the following expression

$$\begin{aligned} M(x)= \left\{ \begin{array}{ll} -\frac{1}{8}(W_0 + 3 W_1) l^2 + \frac{1}{2}(W_0 + W_1) l x -\frac{1}{2} W_0 x^2, &{}\quad 0\le x \le l/2,\\ \left( -\frac{1}{2} l^2 + l x - \frac{1}{2} x^2\right) W_1, &{}\quad l/2<x \le l, \end{array} \right. \end{aligned}$$

(C4)

where $W_0,$ $W_1$ are independent continuous random variables whose p.d.f.’s are known and they are denoted by $f_{W_0}(w_0)$ and $f_{W_1}(w_1) ,$ respectively. We can observe that at $x=l,$ $M(l)=0$.

To obtain the p.d.f. of (C4) and taking into account that it is a piecewise function we need to calculate the p.d.f. separately in the two subdomains.

We fix $x: 0\le x \le l/2$. Now we apply Theorem 1 taking ${\textbf{U}}=\left( W_1, W_0 \right)$ to obtain the p.d.f. of the random vector ${\textbf{V}}=\left( V_1, V_2 \right) ,$ defined by the transformation ${\textbf{r}}:{\mathbb {R}}^2\rightarrow {\mathbb {R}}^2,$ whose components are given by

$$\begin{aligned} v_1= & {} r_1(w_1, w_0) = -\frac{1}{8}(w_0 + 3 w_1) l^2 + \frac{1}{2}(w_0 + w_1) l x -\frac{1}{2} w_0 x^2, \\ v_2= & {} r_2(w_1, w_0) = w_0. \end{aligned}$$

The inverse transformation ${\textbf{s}}: {\mathbb {R}}^{2}\rightarrow {\mathbb {R}}^{2}$ of ${\textbf{r}}$ is given by

$$\begin{aligned} w_1= & {} s_1(v_1,v_2) = \dfrac{v_1 - v_2 \left( -\frac{1}{8} l^2 + \frac{1}{2} l x -\frac{1}{2} x^2 \right) }{\frac{1}{2} l x -\frac{3}{8} l^2}, \\ w_0= & {} s_2(v_1,v_2) = v_2. \end{aligned}$$

Now, we need to calculate the absolute value of the Jacobian of the previous inverse transformation,

$$\begin{aligned} |{\mathcal {J}}| = \left| \frac{\partial s_1(v_1,v_2)}{\partial v_1}\right| = \dfrac{1}{\left| \frac{1}{2} l x -\frac{3}{8} l^2 \right| } \ne 0. \end{aligned}$$

Then, applying the RVT method and taking advantage of the expectation operator we obtain the p.d.f. on the domain $0\le x \le l/2,$

$$\begin{aligned} f_{M(x)}(m)= {\mathbb {E}}_{W_0}\left[ f_{W_1}\left( \dfrac{m-\left( -\frac{1}{8}l^2 + \frac{1}{2} l x - \frac{1}{2} x^2 \right) W_0}{\frac{1}{2} l x - \frac{3}{8} l^2} \right) \dfrac{1}{\left| \frac{1}{2} l x - \frac{3}{8} l^2 \right| } \right] ,\quad 0\le x \le l/2. \end{aligned}$$

Now, we fix $x:l/2<x<l$. In this subdomain, only one random variable appears, so in this case we use the one-dimensional version of the RVT method. Using the same notation as above, we take $U = W_1$ to obtain the p.d.f. of the random variable $V = V_1,$ defined by the transformation $r:{\mathbb {R}}\rightarrow {\mathbb {R}}$ such that

$$\begin{aligned} v_1 = r(w_1) = \left( -\frac{1}{2} l^2 + l x - \frac{1}{2} x^2\right) w_1. \end{aligned}$$

The inverse transformation of r is determined by

$$\begin{aligned} s(v_1):=r^{-1}(v_1) = \dfrac{v_1}{-\frac{1}{2} l^2 + l x - \frac{1}{2} x^2}. \end{aligned}$$

Finally, we need to calculate the derivative of $r^{-1}$ that is given by

$$\begin{aligned} \frac{{{\,\textrm{d}\,}}s}{{{\,\textrm{d}\,}}v_1} = \dfrac{1}{-\frac{1}{2} l^2 + l x - \frac{1}{2} x^2}. \end{aligned}$$

Then applying RVT method, the p.d.f. on $l/2<x<l$ is

$$\begin{aligned} f_{M(x)}(m)= f_{W_1}\left( \dfrac{m}{ - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 } \right) \dfrac{1}{\left| - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 \right| }, \quad l/2<x < l. \end{aligned}$$

In summary, the p.d.f. of the bending moment is given by

$$\begin{aligned} f_{M(x)}(m)= \left\{ \begin{array}{ll} {\mathbb {E}}_{W_0}\left[ f_{W_1}\left( \dfrac{m-\left( -\frac{1}{8}l^2 + \frac{1}{2} l x - \frac{1}{2} x^2 \right) W_0}{\frac{1}{2} l x - \frac{3}{8} l^2} \right) \dfrac{1}{\left| \frac{1}{2} l x - \frac{3}{8} l^2 \right| } \right] , &{}\quad 0\le x \le l/2,\\ f_{W_1}\left( \dfrac{m}{ - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 } \right) \dfrac{1}{\left| - \frac{1}{2} l^2 + l x - \frac{1}{2} x^2 \right| }, &{}\quad l/2<x < l. \end{array} \right. \end{aligned}$$