How accurate are the reliability assessments conducted by the finite element method?

Abstract

The present paper investigates the accuracy of the finite element method (FEM) in stochastic setting. The performance of the FEM for solving the transversal vibration eigenvalue problem of a uniform, homogeneous beam in presence of uncertainties is considered aiming to establish how accurate the method is in predicting the beam’s reliability as well as its probability of failure. An explicit solution is first provided for the approximate fundamental frequency of the beam as a function of the number of elements, for different boundary conditions when the mesh is uniform along the length of the beam allowing an analytical evaluation of the structural reliability and the probability of failure when, e.g., the random uncertainty in the Young modulus of the beam is considered. The exact solution of the vibration problem derived within Bernoulli-Euler beam theory is exploited to evaluate the actual reliability as well as the actual probability of failure which, being compared with required reliability or allowed probability of failure thresholds, permits to verify the accuracy of the FEM in the probabilistic context and to warn about “unreliability of reliability conclusions”.

Introduction

The finite element method (FEM) is almost a universal method to provide approximate solutions to various problems in structural engineering and it has become the dominant engineering simulation tool. According to Oden [1, 2], “Finite elements: Perhaps no other family of approximation methods has had a greater impact on the theory and practice of numerical methods during the twentieth century. Finite element methods have now been used in virtually every area of engineering that can make use of models of nature characterized by partial differential equations”. Its accuracy in deterministic setting has been examined in detail in several investigations, namely, by Pin and Pian [3], Johnson and McLay [4], Babuška [5], and Babuška and Strouboulis [6].

As Mali et al. [7], state: “A mathematical model always represents an “abridged’ version of a physical object, so that the modeling error is always greater than zero”. Moreover, “Approximation errors arise when continual (differential) models are replaced by finite dimensional models”. We deal with the latter error. Additionally, it is not straightforward to extend the considerations on the accuracy of the FEM to the case of scenarios in which the considered problem is affected by uncertain parameters.

Generally, both analysis and design of engineering systems strongly rely on predictions from the numerical models (e.g., the finite element analysis in this context), while the accuracy of the numerical results depends on the closeness between the approximate representation and the real system. However, the feasibility and applicability, as well as the confidence level of a numerical model, such as the FEM, are often challenged by the presence of several uncertainties widely present in engineering practice and deriving from various sources, such as heterogeneity of materials, variability in measurement, and the lack of data.

In the standard FEM, due to its nature as a deterministic analysis tool, all the input variables (based generally on the mean values of the system parameters) are uniquely specified and the output solutions are also uniquely resolved in the form of constant values not representative of all possible scenarios for the system response. Nevertheless, the remarkable influence of the unavoidable intrinsic uncertainties on the behavior of a structural system has motivated the scientific community to recognize the need to resort to stochastic approaches to provide a real prediction of its mechanical response as well as its optimum design.

Issues related to uncertainty quantification, influence of probabilistic model parameters into the reliability of a given structural system and procedures of computational stochastic mechanics have received considerable attention for some time now [8].

In the context of engineering analysis, research efforts on uncertainty quantification are mainly concentrated in two largely research fields [9, 10–11], namely the "Structural Reliability Analysis" and the "Stochastic Finite Element Method" (SFEM).

The structural reliability analysis aims to provide a rational framework for addressing uncertainties in structural analysis so that the design can be more objective. In the stochastic framework, this analysis aims to measure the structural reliability, namely the probability of the structure to perform its mission, or alternatively to evaluate its counterpart, namely its probability of failure.

Despite the simplicity of the definition of the probability of failure, its exact evaluation through direct integration is often very difficult for practical problems, where for real problems the dimension of the integral is usually high and the limit state surface has complicated shape.

On the other hand, in order to account for the various input uncertainties encountered in engineering practice and their effects on the system response, researchers have been trying to extend the standard deterministic FEM into SFEM by incorporating random variables into the mathematical and computational formulations. Specifically, the SFEM is applied for solving stochastic (static and dynamic) problems involving finite elements whose properties are random. In line with the classical FEM, from a mathematical point of view, SFEM can be seen as a powerful tool for the solution of stochastic partial differential equations; this aspect has been treated in several publications with a detailed analysis on convergence and error estimation problems [12, 13, 14, 15, 16–17].

Generally, the first fundamental step in the analysis of uncertain systems (valid for both reliability analysis and SFEM) concerns the representation of the uncertainties (e.g., mechanical and geometric system properties) through the discretization of the stochastic fields [18, 19, 20, 21–22] used to describe such uncertainties and consisting in their approximation (replacement) by a finite number of random variables. The second step of the analysis regards the propagation of uncertainty through the system and the assessment of its stochastic response; this is the most important issue in stochastic mechanics and it is mainly dealt with the proper application of SFEM in a wide variety of engineering problems and in its main variants (see for further details [10] and all the references in it].

In the framework of the analysis of systems in presence of uncertainties, the present paper aims to examine the issue of the accuracy of the finite element method for structures with uncertain parameters from a quite different perspective. The main idea is first to address the problem by applying the standard FEM. Then, by describing an uncertain parameter as a random variable, the evaluation of the structural reliability (as well as the probability of failure) is performed.

Reasonable exploitation of a structure is guaranteed if its reliability is extremely close to unity or alternatively its probability of failure is extremely small, namely, close to zero. The question arises, therefore, about the very ability of the FEM to predict such small probability values with sufficient accuracy. The authors conducted similar analyses on the finite difference method studying the reliability of the reliability calculations in the buckling context [23], longitudinal vibration of a bar [24], and transversal vibration of a beam [25]. Accuracy of the finite difference method in stochastic setting has discussed by Grigoriu and Khater [26]. However, to the best of authors’ knowledge, the FEM in stochastic setting as proposed in this analysis was not considered from the point of view of the convergence.

In the present study, to verify the accuracy of the FEM, first the natural frequency of a uniform beam (with a specific attention on a simply supported beam) known exactly by the Euler–Bernoulli beam theory, has been approximately calculated by properly applying the FEM returning an analytical expression in terms of the beam parameters and the number of segments adopted for the mesh. By introducing the probabilistic characterization of, e.g., the Young modulus parameter (modeled for simplicity as a continuous random variable with assigned probability density function) and taking full advantage of the derived explicit solution, the expressions for the so-called actual reliability and the actual probability of failure are derived. That is, noting that generally a structure performs satisfactorily if its reliability in not less than a codified reliability value, or-alternatively- if the probability of failure is less than an its tolerable level, a design problem can be solved by selecting some design parameter (e.g., the length of the beam) depending on the number of mesh elements and the value of the reliability (or probability of failure) threshold. The actual values of reliability or probability of failure are lastly calculated considering their exact expressions where the design parameter is substituted by the value provided via the FEM approximate analysis.

The results in terms of percentage error in reliability and probability of failure evaluated with respect to the discretization parameter for various values of codified reliabilities and codified levels of probability of failure allow to draw some unexpected conclusions on the accuracy of FEM.

FEM analysis via Seide’s solution

In this section an approach via FEM applied to the problem of a column buckling as discussed by Seide [27] has been developed to analyze the transversal vibrations problem of a beam. Thomée [28] discussed both finite difference and finite element method in terms of “simple model situations”. We follow his methodological advice. Specifically, we consider the transverse vibrations of a uniform and homogeneous Euler–Bernoulli beam governed by the following differential equation:

1

$$EI\frac{{\partial }^{4}w\left(x,t\right)}{\partial x^4}+\rho A\frac{{\partial }^{2}w\left(x,t\right)}{\partial t^2}=0$$where E is the Young modulus of elasticity, I is the inertia moment, A is the cross-sectional area, $\rho$ is the mass density, x is the axial coordinate, t is the time and w the transverse displacement.

For beam in free vibrations, by setting $w\left(x,t\right)=W\left(x\right)sin\omega t$, Eq. (1) can be rewritten as:

2

$$EI\frac{{\text{d}}^{4}W\left(x\right)}{\text{d}{x}^4}-\rho A{\omega }^{2}W\left(x\right)=0$$

Solution of Eq. (2) can be provided by resorting to FEM based on the principle of minimum potential energy. The variational principle for the beam in transversal vibrations leads to:

3

$$\delta \underset{0}{\overset{L}{\int }}\frac{1}{2}\left[EI{\left(\frac{{\text{d}}^{2}W}{\text{d}{x}^2}\right)}^2-\rho A{\omega }^2{W}^2\right]\text{d}x=0$$ where in the notation the dependency of W by x is omitted for simplicity.

By dividing the length of the beam into N elements, the deflection function of each element may be expressed as a cubic function whose coefficients in its expression are given in terms of nodal transversal displacements and slopes at each end of the segment. Thus, by introducing the local axial coordinate $\xi =x\left./\phantom{x{\ell }_i}\right){\ell }_i$ being ${\ell }_i$ the length of the i-th element, the deflection function takes the following expression:

4

$$W=W_{i-1}\left(1-3{\xi }^2+2{\xi }^3\right)+W_i{\xi }^2\left(3-2\xi \right)+W_{i-1}\prime {\ell }_i\xi {\left(1-\xi \right)}^2-W_i\prime {\ell }_i{\xi }^2\left(1-\xi \right);0<\xi <1$$where $W_{i-1}$ and $W_{i-1}\prime$ represent deflection and slope at $\xi =0$ while $W_i$ and $W_i\prime$ represent deflection and slope at $\xi =1$, respectively.

Equation (3) is now rewritten as:

5

$$\delta \sum _{i=0}^N\underset{0}{\overset{1}{\int }}\frac{1}{2}\left[{\left(\frac{{\text{d}}^{2}W}{\text{d}{\xi }^2}\right)}^2-\frac{\rho A{\omega }^2{\ell }_i^4}{\mathit{EI}}{W}^2\right]\text{d}\xi =0$$

By performing the variation in Eq. (5) for all the segments with same length ${\ell }_i=\ell$ and by considering the deflection expression in Eq. (4), the two following equations are provided:

6a

$$\begin{array}{c}-\left(W_{i-1}+W_{i+1}\right)\left(12+\frac{9}{70}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)+2W_i\left(12-\frac{13}{35}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)+\\ +\left(W_{i+1}\prime \ell -W_{i-1}\prime \ell \right)\left(6+\frac{13}{420}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)=0\\ \end{array}$$

6b

$$\begin{array}{c}\left(W_{i-1}-W_{i+1}\right)\left(6+\frac{13}{420}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)+2W_i\prime \ell \left(4-\frac{1}{150}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)+\\ +\left(W_{i-1}\prime \ell +W_{i+1}\prime \ell \right)\left(2+\frac{1}{140}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)=0\\ \end{array}$$

The system in Eqs. (6a) and (6b) is solved by expressing $W_i$ and $W_i\prime$ in terms of $\lambda$ as

7a

$$W_i=B{\lambda }^i$$

7b

$$W_i^{\prime }=C{\lambda }^i$$

By substituting Eqs. (7a) and (7b) in Eqs. (6a) and (6b) the nontrivial solution leads to the equation in $\lambda$ expressed by:

8

$$a{\lambda }^4-4b{\lambda }^3+2c{\lambda }^2-4b\lambda +a=0$$where the following positions hold:

9

$$a=1+\frac{1}{420}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2+\frac{1}{302400}{\left(\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)}^2;$$

10

$$b=1+\frac{37}{840}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2+\frac{1}{16800}{\left(\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)}^2;$$

11

$$c=3-\frac{137}{420}\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2+\frac{131}{302400}{\left(\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2\right)}^2$$

After some manipulations, solution of Eq. (8) is provided by:

12a, b

$$\begin{array}{c}{\lambda }_{1,2}=\frac{b}{a}-\frac{1}{2}\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\mp \sqrt{\frac{1}{4}{\left(\frac{2b}{a}-\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\right)}^2-1}\\ {\lambda }_{3,4}=\frac{b}{a}+\frac{1}{2}\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\mp \sqrt{\frac{1}{4}{\left(\frac{2b}{a}+\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\right)}^2-1}\\ \end{array}$$where a, b and c take the expressions given in Eqs. (9), (10) and (11), respectively.

By defining

13

$$cos\vartheta =\frac{1}{2}\left(2\frac{b}{a}-\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\right)$$the general solution for $W_i$ can be written as follows:

14

$$W_i=D_{1}cosi\vartheta +D_{2}sini\vartheta +D_3{\lambda }_3^i+D_4{\lambda }_4^i$$with $D_1,D_2,D_3$ and $D_4$ arbitrary constants of integration to determine by fixing the boundary conditions for the beam.

By analyzing i.e.a simply supported beam the boundary conditions read:

15

$$W_0=W_N=0;W_{-1}=-W_1;W_{N+1}=-W_{N-1}$$and taking into account Eq. (14) the following four equations are provided:

16a

$$D_1+D_3+D_4=0$$

16b

$$D_{1}cosN\vartheta +D_{2}sinN\vartheta +D_3{\lambda }_3^N+D_4{\lambda }_4^N=0$$

16c

$$2D_{1}cos\vartheta +D_3\left(\frac{1}{\lambda _3^{}}+\lambda _3^{}\right)+D_4\left(\frac{1}{\lambda _4^{}}+\lambda _4^{}\right)=0$$

16d

$$2D_{1}cos\vartheta cosN\vartheta +2D_{2}cos\vartheta sinN\vartheta +D_3{\lambda }_3^N\left(\frac{1}{\lambda _3^{}}+\lambda _3^{}\right)+D_4{\lambda }_4^N\left(\frac{1}{\lambda _4^{}}+\lambda _4^{}\right)=0$$

The nontrivial solution requires the determinant of the coefficients equal to zero that is:

17

$$\left({\lambda }_3^N-{\lambda }_4^N\right)\left(1+{\lambda }_3^2-2\lambda _3^{}cos\vartheta \right)\left(1+{\lambda }_4^2-2\lambda _4^{}cos\vartheta \right)sinN\vartheta =0$$

Bearing in mind Eq. (12b) and by introducing the positions

18a, b

$$\text{d}=\frac{b}{a}\text{and}f=\frac{1}{2}\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}$$

Equation (17) can be rewritten as

19

$$4\left[{\left(d+f-\sqrt{{(d+f)}^2-1}\right)}^N-{\left(d+f+\sqrt{{(d+f)}^2-1}\right)}^N\right]{\left(d+f-cos\vartheta \right)}^{2}sinN\vartheta =0$$

The satisfaction of Eq. (19) involves:

20

$${\left(d+f-\sqrt{{(d+f)}^2-1}\right)}^N={\left(d+f+\sqrt{{(d+f)}^2-1}\right)}^N$$or

21

$${\left(d+f-cos\vartheta \right)}^2=0$$or

22

$$sinN\vartheta =0$$

Equality in Eq. (20) implies

23

$${(d+f)}^2-1={\left(\frac{b}{a}+\frac{1}{2}\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\right)}^2-1=0$$whose solution, keeping in mind Eqs. (9),(10) and (11) returns

24

$$\frac{\rho A}{\mathit{EI}}{\ell }^4{\omega }^2=0$$that is $\omega =0$, a meaningless result to be discarded.

Equation (21) requires

25

$$cos\vartheta =d+f=\frac{b}{a}+\frac{1}{2}\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}$$but the right side takes values greater than 1 and it is in contrast with the position in Eq. (13) and therefore this result must also be discarded.

Then the condition of determinant of the coefficients equal to zero reduces to Eq. (22) which returns:

26

$$N\vartheta =k\pi k=1,2,3,...$$

Evaluation of $cos\vartheta$ leads to

27

$$cos\frac{k\pi }{N}=\frac{1}{2}\left(2\frac{b}{a}-\sqrt{2+{\left(\frac{2b}{a}\right)}^2-\frac{2c}{a}}\right)$$

By substituting in Eq. (27) the expressions for a,b and c (see Eqs. 9, 10, 11), setting k = 1 for the first natural frequency and bearing in mind that the length $\ell$ of each segment is equal to $L\left./\phantom{LN}\right)N$ the beam fundamental frequency can be calculated as

28

$${\omega }_1^2=\frac{{\pi }^4}{L^4}\frac{\mathit{EI}}{\rho A}g(N){\left(\frac{N}{\pi }\right)}^4$$with

29

$$g(N)=\frac{120\left(-411-222\text{cos}\frac{\pi }{N}+3\text{Cos}\frac{2\pi }{N}+\sqrt{6}\sqrt{30380+33026\text{cos}\frac{\pi }{N}+2723\text{cos}\frac{2\pi }{N}+22\text{cos}\frac{3\pi }{N}-\text{cos}\frac{4\pi }{N}}\right)}{72\text{cos}\frac{\pi }{N}-\text{cos}\frac{2\pi }{N}-131}$$

Equation (28) provides a solution for the fundamental frequency for the analyzed beam with an explicit dependence on the mechanical and geometrical beam parameters as well as on the number of the beam elements employed in the FEM. The latter explicit solution allows performing the analysis conducted in the next section.

It is worth to note that the limit of $g(N){\left(N\left./\phantom{N\pi }\right)\pi \right)}^4$ when N goes to infinity is equal 1 and the exact expression of the frequency for a simply supported beam in free vibrations is provided.

Probabilistic analysis for random elastic modulus

Let us examine the case where the elastic Young modulus E of the beam is described as a continuous random variable with assigned probability density function (PDF). All the other parameters involved in the analysis are treated as deterministic quantities.

The analysis is aimed to avoid resonance phenomenon that is the natural frequency of the beam must be less than an excitation frequency value $\omega _0^{}:$

30

$$\omega _1^{}<\omega _0^{}$$

Since the elastic Young modulus E is treated as a random variable, it follows that the fundamental frequency $\omega _1^{}$ as reported in Eq. (28) is in turn a random variable too and labeled in the following by the capital letter $\mathrm{\Omega }_1^{}$ as generally agreed with the random variables notation.

To analyze the possible values of E satisfying the inequality in Eq. (30), the reliability R is introduced as

31

$$R\text{= Prob}\left(\mathrm{\Omega }_1^{}<\omega _0^{}\right)$$that is the probability of the event reported in Eq. (30).

Starting from its definition, in Eq. (31) the reliability R can be written in terms of the approximate fundamental frequency value reported in Eq. (28) leading to

32

$$R_{\text{approx}}\text{= Prob}\left(\frac{{\pi }^4}{L^4}\frac{\mathit{EI}}{\rho A}g(N){\left(\frac{N}{\pi }\right)}^4<{\omega }_0^2\right)$$

Under the previous assumption that all the involved parameters are considered deterministic quantities with the exception of the Young modulus, the approximate reliability can be expressed as

33a

$$R_{\text{approx}}\text{= Prob}\left(E<\frac{L^4}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\frac{1}{g(N)}{\left(\frac{\pi }{N}\right)}^4\right)$$or alternatively by introducing the probability distribution function of E,$F_E\left(e\right)$, as:

33b

$$R_{\text{approx}}=F_E\left(\frac{L^4}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\frac{1}{g(N)}{\left(\frac{\pi }{N}\right)}^4\right)$$

Based on the consideration that a structure performs satisfactorily if the reliability is greater or equal than a codified reliability value $r_0$, the design problem can be expressed as:

34

$$R\ge r_0,0<r_0\le 1$$

Alternatively, one can introduce the probability of failure $P_f$:

35

$$P_f=1-R$$and define the design problem as

36

$$P_f\le p_0$$where $p_0$ represents the tolerable level of probability of failure.

As highlighted by Eqs. (34–36), the knowledge of the reliability R allows to solve the design problem for the beam. Moreover, for the case under exam, the exact expression of the natural frequency for the simply supported beam is given by the well-known expression

37

$${\omega }_{\text{exact}}=\frac{{\pi }^2}{L^2}\sqrt{\frac{\mathit{EI}}{\rho A}}$$which allows to calculate the exact reliability $R_{\text{exact}}$ for the beam.

Let us assume, e.g., the length of the beam L as a design parameter: its expression can be derived by the approximate analysis once the PDF for the random variable E is assigned. It is worth to note that the design parameter $L_{\text{approx}}$ depends on the number of segment N and the specified value for $r_0$ that is $L_{\text{approx}}=L\left(N,r_0\right).$

The substitution of the approximate parameter $L_{\text{approx}}$ into the reliability expression deduced by considering ${\omega }_{\text{exact}}$ in Eq. (37) provides a general expression for the so-called actual reliability $R_{\text{actual}}$ according to N and $r_0$ values. The comparison between the actual reliability and the required reliability $r_0$ allows to verify the accuracy of the FEM in stochastic setting.

To clarify the aim of the described analysis two numerical examples are reported in the following.

Results for exponentially distributed beam elastic modulus

Let us assume that the random variable E is probabilistically characterized by, e.g., an exponential distribution given as

38

$$f_E\left(e\right)=\left\{\begin{array}{c}0e<0\\ aexp\left[-ae\right]e\ge 0,a>0\\ \end{array}\right)$$with mean value $<E>=1\left./\phantom{1a}\right)a$ and variance ${\sigma }_E^2=1\left./\phantom{1a^2}\right)a^2$, where the symbol $<·>$ stands for mathematical expectation.

The expression of the approximate reliability $R_{\text{approx}}^{\text{(exp)}}$ for the exponential distribution of E, reads:

39

$$R_{\text{approx}}^{\text{(exp)}}=1-exp\left[-\frac{1}{<E>}\frac{L^4}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\frac{1}{g(N)}{\left(\frac{\pi }{N}\right)}^4\right]$$

In design setting, by posing $R_{\text{approx}}^{\text{(exp)}}=r_0$ and by performing some manipulations, the design parameter $L_{\text{approx}}^{\text{(exp)}}=L\left(N,r_0\right)$ is provided as:

40

$$L_{\text{approx}}^{\text{(exp)}}=\pi \sqrt[4]{<E>\frac{I}{\rho A{\omega }_0^2}g(N)ln\left(\frac{1}{1-r_0}\right)}\left(\frac{N}{\pi }\right)$$

Bearing in mind Eq. (37), the exact reliability can be expressed as

41

$$R_{\text{exact}}=\text{Prob}\left(\frac{{\pi }^4}{L^4}\frac{\mathit{EI}}{\rho A}<{\omega }_0^2\right)$$or, based on Eq. (38), as:

42

$$R_{\text{exact}}=1-exp\left[-\frac{1}{<E>}\frac{L^4}{{\pi }^4}\frac{{\omega }_0^2\rho A}{I}\right]$$

Finally, the actual reliability $R_{\text{actual}}^{\text{(exp)}}$ can be achieved by considering the approximate parameter in the $L_{\text{approx}}^{\text{(exp)}}=L\left(N,r_0\right)$ expression of $R_{\text{exact}}$ as:

43

$$R_{\text{actual}}^{\text{(exp)}}=R_{\text{actual}}^{\text{(exp)}}\left(N,r_0\right)={\left(R_{\text{exact}}\right|}_{L=L_{\text{approx}}^{(exp)}}=1-exp\left[-\frac{1}{<E>}\frac{L_{\text{actual}}^{{\text{(exp)}}^4}}{{\pi }^4}\frac{{\omega }_0^2\rho A}{I}\right]$$

The substitution of the expression of $L_{\text{approx}}^{\text{(exp)}}=L\left(N,r_0\right)$ (see Eq. 40) in Eq. (43) leads to

44

$$R_{\text{actual}}^{\text{(exp)}}=R_{\text{actual}}^{\text{(exp)}}\left(N,r_0\right)=1-exp\left[g(N)\left(\frac{N^4}{{\pi }^4}\right)ln(1-r_0^{})\right]=1-(1-r_0^{})_{}^{g(N){\left(N\left./\phantom{N\pi }\right)\pi \right)}^4}$$

The evaluation of $R_{\text{actual}}^{\text{(exp)}}$ via Eq. (44) returns values higher than the codified value $r_0$. The following graphics (Figs. 1a–d) show the trend of the percentage error ${\epsilon }^{(exp)}\left(\%\right)$ defined as:

45

$${\epsilon }^{(exp)}\left(\%\right)=\frac{\left|r_0^{}-R_{\text{actual}}^{\text{(exp)}}\right|}{r_0}\times 100$$

[See PDF for image]

Fig. 1

Percentage error ${\epsilon }^{(exp)}\left(\%\right)$ in reliability evaluation versus the discretization parameter N for various codified reliabilities $r_0$: a$r_0=0.90$; b$r_0=0.99$; c$r_0=0.999$ and d$r_0=0.9999$

For $r_0=0.90$ and N = 1 the error is ${\epsilon }^{(exp)}\left(\%\right)=0.37\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.903288$) decreasing to ${\epsilon }^{(exp)}\left(\%\right)=0.0055\%$$(R_{\text{actual}}^{\text{(exp)}}=0.900049)$ for N = 5 and ${\epsilon }^{(exp)}\left(\%\right)=0.0003\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.900003$) for N = 10. By selecting $r_0=0.99$ the error goes to ${\epsilon }^{(exp)}\left(\%\right)=0.065\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.990647$) for N = 1 to ${\epsilon }^{(exp)}\left(\%\right)=0.00099\%$$\left(R_{\text{actual}}^{\text{(exp)}}=0.99001\right)$ for N = 5. For $r_0=0.999$ and $r_0=0.9999$ the evaluated errors are ${\epsilon }^{(exp)}\left(\%\right)=0.0096\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.999095$) and ${\epsilon }^{(exp)}\left(\%\right)=0.00125\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.999913$) for N = 1 and ${\epsilon }^{(exp)}\left(\%\right)=0.00015\%$$\left(R_{\text{actual}}^{\text{(exp)}}=0.999001\right)$ and ${\epsilon }^{(exp)}\left(\%\right)=0.00002\%$ ($R_{\text{actual}}^{\text{(exp)}}=0.9999$) for N = 5, respectively. Analyzing the latter results as well as by the inspection of the graphics it follows that for a fixed value of N, the percentage error decreases when the codified value of reliability $r_0$ increases.

Starting from its definition (see Eq. 35), the actual probability of failure can be expressed as:

46

$$P_{f\text{,actual}}^{(exp)}=p_0^{\gamma }$$with

47

$$\gamma =g(N)(N\left./\phantom{N{\pi )}^4}\right){\pi )}^4$$where the expression of g(N) is reported in Eq. (29). It is worth to note that $\gamma$ in Eq. (47) takes values slightly higher than unity and very close to unity with the increasing number of segments N: that is $\gamma =1.01452$ for N = 1, $\gamma =1.00021$ for N = 5 and $\gamma =1.00001$ for N = 10.

Graphics in Fig. 2a–d show the variation of the actual probability of failure $P_{f\text{,actual}}^{(exp)}$ versus the number of segments N for fixed levels of the codified value of probability of failure $p_0$; the actual probability of failure is always smaller than the requested value. It is worth to note that different results have been previously obtained by solving the same problem via finite difference method (see [25]) where the curves representing the actual probabilities of failure for the different codified thresholds are always in the unsafe side.

[See PDF for image]

Fig. 2

Actual probability of failure $P_{f\text{,actual}}^{(exp)}$ versus the discretization parameter N or various codified levels $p_0$: a$p_0=0.1$; b$p_0=0.01$; c$p_0=0.001$ and d$p_0=0.0001$

For a fixed value $p_0=0.1$ the values of the probability of failure obtained are $P_{f\text{,actual}}^{(exp)}=0.967118p_0$ for N = 1 and $P_{f\text{,actual}}^{(exp)}=0.999507p_0$ for N = 5. For $p_0=0.01$ one gets $P_{f\text{,actual}}^{(exp)}=0.935317p_0$ and $P_{f\text{,actual}}^{(exp)}=0.999015p_0$ for N = 1 and N = 5, respectively. The values $P_{f\text{,actual}}^{(exp)}=0.904562p_0$ (N = 1) and $P_{f\text{,actual}}^{(exp)}=0.998522p_0$ (N = 5) are provided for $p_0=0.001$ while $P_{f\text{,actual}}^{(exp)}=0.874818p_0$ (N = 1) and $P_{f\text{,actual}}^{(exp)}=0.99803p_0$ (N = 5) are provided for $p_0=0.0001.$ For a fixed value of $p_0$, with the increasing number of N, the value of $P_{f\text{,actual}}^{(exp)}$ increases, reaching the $p_0$ value at a low number of N (namely N = 6/7).

Results for normal distributed beam elastic modulus

Let us analyze the case in which the random variable E is characterized from a probabilistic point of view by a normal or Gaussian distribution; the probability density function is thus expressed as:

48

$$f_E\left(e\right)=\frac{1}{\sqrt{2\pi }{\sigma }_E}exp\left[-\frac{e^2}{2{\sigma }_E^2}\right]$$with zero mean value and variance ${\sigma }_E^2$.

For the case under exam the approximate reliability $R_{\text{approx}}^{\text{(norm)}}$ is provided by

49

$$R_{\text{approx}}^{\text{(norm)}}=\frac{1}{2}\text{erfc}\left[-\frac{1}{\sqrt{2}{\sigma }_E}\frac{L^4}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\frac{1}{g(N)}{\left(\frac{\pi }{N}\right)}^4\right]$$where erfc[z] gives the complementary error function as

50

$$\text{erfc}\left[z\right]=1-\text{erf}\left[z\right]=1-\frac{2}{\sqrt{\pi }}\underset{0}{\overset{z}{\int }}exp\left[-t^2\right]dt$$with erf[z] in turn denoting the error function namely the integral of the Gaussian distribution.

By following the steps reported in the previous section, the design parameter $L_{\text{approx}}^{\text{(norm)}}=L\left(N,r_0\right)$ is provided by assuming $R_{\text{approx}}^{\text{(norm)}}=r_0$ in Eq. (49) and solving for L:

51

$$L_{\text{approx}}^{\text{(norm)}}=\pi \sqrt[4]{-\sqrt{2}{\sigma }_E\frac{I}{\rho A{\omega }_0^2}g(N)\text{erfcinv}\left[2r_0\right]}\left(\frac{N}{\pi }\right)$$where in Eq. (51) $\text{erfcinv}\left[z\right]$ denotes the inverse complementary error function obtained as the solution for z in $s=\text{erfc}\left[z\right]$.

Recalling Eq. (37) and Eq. (41), due to the probability distribution on Eq. (48), the exact reliability can be rewritten as

52

$$R_{\text{exact}}=\frac{1}{2}\text{erfc}\left[-\frac{1}{\sqrt{2}{\sigma }_E}\frac{L^4}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\right]$$and accordingly, the actual reliability $R_{\text{actual}}^{\text{(norm)}}$ is obtained by substituting the approximate design parameter $L_{\text{approx}}^{\text{(norm)}}=L\left(N,r_0\right)$ in Eq. (52) for the exact reliability as:

53

$$R_{\text{actual}}^{(\text{norm})}=R_{\text{actual}}^{(\text{norm})}\left(N,r_0\right)={\left(R_{\text{exact}}\right|}_{L=L_{\text{approx}}^{\text{(norm)}}}=\frac{1}{2}\text{erfc}\left[-\frac{1}{\sqrt{2}{\sigma }_E}\frac{L_{\text{approx}}^{{\text{(norm)}}^4}}{{\pi }^4}\frac{\rho A{\omega }_0^2}{I}\right]$$

Bearing in mind Eq. (51), Eq. (53) can be rewritten as

54

$$R_{\text{actual}}^{(\text{norm})}=R_{\text{actual}}^{(\text{norm})}\left(N,r_0\right)={\left(R_{\text{exact}}\right|}_{L=L_{\text{approx}}^{(\text{norm})}}=\frac{1}{2}\text{erfc}\left[g(N)\left(\frac{N^4}{{\pi }^4}\right)\text{erfcinv}[2r_0]\right]$$

As for the case of $R_{\text{actual}}^{(\text{exp})}$, calculations of $R_{\text{actual}}^{(\text{norm})}$ through Eq. (54) return values higher than the codified value $r_0$. The percentage error ${\epsilon }^{\text{(norm)}}\left(\%\right)$ for the case of normal distribution for the random Young modulus E variable is defined as

55

$${\epsilon }^{\text{(norm)}}\left(\%\right)=\frac{\left|r_0^{}-R_{\text{actual}}^{(\text{norm})}\right|}{r_0}\times 100$$and its trend is depicted in the graphics reported in Fig. 3a–d for increasing value of $r_0$.

[See PDF for image]

Fig. 3

Percentage error ${\epsilon }^{\text{(norm)}}\left(\%\right)$ in reliability evaluation versus the discretization parameter N for various codified reliabilities $r_0$: a$r_0=0.90$; b$r_0=0.99$; c$r_0=0.999$ and d$r_0=0.9999$

In this second case study and for $r_0=0.90$ the error is ${\epsilon }^{\text{(norm)}}=0.36\%$($R_{\text{actual}}^{(\text{norm})}=0.903227$) for N = 1 decreasing to ${\epsilon }^{\text{(norm)}}=0.0053\%$$(R_{\text{actual}}^{(\text{norm})}=0.900048)$ for N = 5 and ${\epsilon }^{\text{(norm)}}=0.00033\%$$\left(R_{\text{actual}}^{(\text{norm})}=0.900005\right)$ for N = 10. Assuming $r_0=0.99,$ the error goes to ${\epsilon }^{\text{(norm)}}=0.087\%$($(R_{\text{actual}}^{(\text{norm})}=0.990866)$ for N = 1 (an error ${\epsilon }^{\text{(exp)}}=0.065\%$ for the previous case of exponential distribution) to ${\epsilon }^{\text{(norm)}}=0.000648\%$$\left(R_{\text{actual}}^{(\text{norm})}=0.990013\right)$ for N = 5.

By considering $r_0=0.999$ and $r_0=0.9999$ the errors calculated for N = 1 are ${\epsilon }^{\text{(norm)}}=0.014\%$$(R_{\text{actual}}^{(\text{norm})}=0.999141)$ and ${\epsilon }^{\text{(norm)}}=0.0019\%$ ($R_{\text{actual}}^{(\text{norm})}=0.999919$) and ${\epsilon }^{\text{(norm)}}=0.00022\%$$\left(R_{\text{actual}}^{(\text{norm})}=0.999002\right)$ and ${\epsilon }^{\text{(norm)}}=0.00003\%$ ($R_{\text{actual}}^{(\text{norm})}==0.9999$) for N = 5, respectively.

The analysis of the previously reported results as well as the trend of the graphics lead again that for a fixed value of N, the percentage error decreases when the codified value of reliability $r_0$ increases.

Based on its definition reported in Eq. (35), the actual probability of failure $P_{f\text{,actual}}^{(\text{norm})}$ can be expressed as:

56

$$P_{f\text{,actual}}^{(\text{norm})}=1-R_{\text{actual}}^{(\text{norm})}=1-\frac{1}{2}\text{erfc}\left[g(N)\left(\frac{N^4}{{\pi }^4}\right)\text{erfcinv}[2(1-p_0)]\right]$$

The trend of the actual probability of failure $P_{f\text{,actual}}^{(\text{norm})}$ versus the number of segments N for fixed levels of the codified value of probability of failure $p_0$ are reported in Fig. 4a–d; also for this second case the actual probability of failure is always smaller than the requested value.

[See PDF for image]

Fig. 4

Actual probability of failure $P_{f\text{,actual}}^{(\text{norm})}$ versus the discretization parameter N or various codified levels $p_0$: a$p_0=0.1$; b$p_0=0.01$; c$p_0=0.001$ and d$p_0=0.0001$

By fixing the value $p_0=0.1$ the values of the probability of failure obtained are $P_{f\text{,actual}}^{(\text{norm})}=0.096773$ for N = 1 and $P_{f\text{,actual}}^{(\text{norm})}=0.0999519$ for N = 5. For $p_0=0.01$ the two values $P_{f\text{,actual}}^{(\text{norm})}=0.00913432$ and $P_{f\text{,actual}}^{(\text{norm})}=0.00998674$ are obtained for N = 1 and N = 5, respectively. The values $P_{f\text{,actual}}^{(\text{norm})}=0.0008589$ (N = 1) and $P_{f\text{,actual}}^{(\text{norm})}=0.0009978$ (N = 5) are provided for $p_0=0.001$ while $P_{f\text{,actual}}^{(\text{norm})}=0.00008064$ (N = 1) and $P_{f\text{,actual}}^{(\text{norm})}=0.00009968$ (N = 5) are provided for $p_0=0.0001.$

The choice of the probability distribution for the random elastic modulus of the beam is irrelevant on the results: analogous qualitative results are obtained by selecting, e.g., a Rayleigh distribution as well as a uniform distribution for the random selected parameter.

Conclusions

The paper investigates the accuracy of the FEM in the transversal vibrations of a beam taking full advantage of a derived closed-form expression in terms of the discretization parameter as well as the parameters characterizing the beam. A probabilistic analysis is conducted by taking into account the uncertainty in the material property (namely the elastic modulus of the beam) and drawing some interesting considerations observing the results in terms of reliability as well as probability of failure in the design beam problem. It is worth highlighting the unexpected results in terms of reliability calculation where for a fixed number of elements of the mesh, the error decreases with an increase of the codified value of reliability. This means that the behavior of the structure (the beam in the case under exam) improves for higher values of reliability set in the design process providing a remarkable result. The probability of failure converges to its codified threshold increasing the number of the mesh elements maintaining its value always smaller than it: the actual value obtained is in the safe side for the design of the beam and it represents a good new in the evaluation of the accuracy of FEM.

Acknowledgements

This paper is dedicated to the blessed memory of Professor John Tinsley Oden.

Author contributions

R.S: Conceptualization, Methodology, Writing—editing. I.E: Conceptualization, Methodology. All authors reviewed the manuscript.

Data availability

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Declarations

Conflict of interest

The authors declare no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Oden, J.T.: Some historic comments on finite elements, In: Proceedings of the ACM conference on History of Scientific and Numeric Computation, pp. 125–130 (1987)

2. Oden J.T.: Historical comments on finite elements. In: A History of Scientific Computing, pp. 152–166 (1990)

3. Pin, T; Pian, THH. The convergence of finite element method in solving linear elastic problems. Int. J. Solids Struct.; 1967; 3, 5 pp. 865-879. [DOI: https://dx.doi.org/10.1016/0020-7683(67)90059-5] 0149.42802

4. Johnson, MW, Jr; McLay, RW. Convergence of the Finite Element Method in the Theory of Elasticity. J. Appl. Mech.; 1968; 35, 2 pp. 274-278. [DOI: https://dx.doi.org/10.1115/1.3601191] 0164.26101

5. Babuška, I. Error bounds for the finite element method. Numer. Math.; 1971; 16, pp. 322-333.288971 [DOI: https://dx.doi.org/10.1007/BF02165003] 0214.42001

6. Babuška, I; Strouboulis, T. The Finite Element Method and Its Reliability; 2001; Oxford, Oxford University Press: [DOI: https://dx.doi.org/10.1093/oso/9780198502760.001.0001] 0995.65501

7. Mali, O; Neittaanmäki, P; Repin, S. Accuracy Verification Methods: Theory and Algorithms; 2014; Dordrecht, Springer: [DOI: https://dx.doi.org/10.1007/978-94-007-7581-7] 1364.65106

8. Oden, JT; Belytschko, T; Babuška, I; Hughes, TJR. Research directions in computational mechanics. Comput. Methods Appl. Mech. Eng.; 2003; 192, pp. 913-922.1954214 [DOI: https://dx.doi.org/10.1016/S0045-7825(02)00616-3] 1025.74503

9. Der Kiureghian, A; KeJ, B. The stochastic finite element method in structural reliability Probabilist. Eng. Mech.; 1988; 3, pp. 83-91.0765.60085

10. Stefanou, G. The stochastic finite element method: past, present and future. Comput. Methods Appl. Mech. Eng.; 2009; 198, 9–12 pp. 1031-1051. [DOI: https://dx.doi.org/10.1016/j.cma.2008.11.007] 1229.74140

11. Aldosary, M; Wang, J; Li, C. Structural reliability and stochastic finite element methods: state-of-the-art review and evidence-based comparison. Eng. Comput.; 2018; 35, 6 pp. 2165-2214. [DOI: https://dx.doi.org/10.1108/EC-04-2018-0157] 1078.81042

12. Babuška, I; Chatzipantelidis, P. On solving elliptic stochastic partial differential equationsComput. Methods Appl. Mech. Eng.; 2002; 191, pp. 4093-4122. [DOI: https://dx.doi.org/10.1016/S0045-7825(02)00354-7] 1019.65010

13. Babuška, I; Tempone, R; ZourarisG, E. Galerkin finite element approximations ofstochastic elliptic partial differential equationsSIAM. J. Numer. Anal.; 2004; 42, pp. 800-825. [DOI: https://dx.doi.org/10.1137/S0036142902418680] 1080.65003

14. Babuška, I; Tempone, R; ZourarisG, E. Solving elliptic boundary value problems with uncertain coefficients by the finite element method: the stochastic formulation. Comput. Methods Appl. Mech. Eng.; 2005; 194, pp. 1251-1294.2121215 [DOI: https://dx.doi.org/10.1016/j.cma.2004.02.026] 1087.65004

15. DebM, K; BabuškaI, M; Tinsley, OJ. Solution of stochastic partial differential equations using Galerkin finite element techniques. Comput. Methods Appl. Mech. Eng.; 2001; 190, pp. 6359-6372.1870425 [DOI: https://dx.doi.org/10.1016/S0045-7825(01)00237-7] 1075.65006

16. Doostan, A; Ghanem, RG; Red-Horse, J. Stochastic model reduction for chaos representations. Comput. Methods Appl. Mech. Eng.; 2007; 196, pp. 3951-3966.2340015 [DOI: https://dx.doi.org/10.1016/j.cma.2006.10.047] 1173.74411

17. Frauenfelder, P; Schwab, C; TodorR, A. Finite elements for elliptic problems with stochastic coefficients. Comput. Methods Appl. Mech. Eng.; 2005; 194, pp. 205-228.2105161 [DOI: https://dx.doi.org/10.1016/j.cma.2004.04.008] 1143.65392

18. Nouy, A; Clément, A; Schoefs, F; Moës, N. An extended stochastic finite element method for solving stochastic partial differential equations on random domains. Comput. Methods Appl. Mech. Eng.; 2008; 197, pp. 4663-4682.2464510 [DOI: https://dx.doi.org/10.1016/j.cma.2008.06.010] 1194.74457

19. Shinozuka, M; Deodatis, G. Simulation of multi-dimensional Gaussian stochastic fields by spectral representation. Appl. Mech. Rev. (ASME); 1996; 49, pp. 29-53. [DOI: https://dx.doi.org/10.1115/1.3101883] 0976.74559

20. Ghanem, R; Spanos, PD. Stochastic Finite Elements: A Spectral Approach; 2003; 2 Berlin, Springer-Verlag: 0722.73080

21. Kleiber, M; Hien, TD. The Stochastic Finite Element Method; 1992; New York, John Wiley & Sons: 0902.73004

22. Papadrakakis, M; Papadopoulos, V. Robust and efficient methods for stochastic finite element analysis using Monte Carlo simulation. Comput. Methods Appl. Mech. Eng.; 1996; 134, pp. 325-340.1412009 [DOI: https://dx.doi.org/10.1016/0045-7825(95)00978-7] 0891.73079

23. Elishakoff, I., Santoro, R.: Reliability of structural reliability estimation. In: Muhanna, R.L., Mullen, R.L. (eds.) Proceedings of the NSF Workshop on Reliable Engineering Computing, pp. 53–64, Georgia University Press (2006)

24. Elishakoff, I; Santoro, R. Error in finite difference based probabilistic dynamic analysis: analytical evaluation. J. Sound Vib.; 2005; 281, 3–5 pp. 1195-1206. [DOI: https://dx.doi.org/10.1016/j.jsv.2004.03.066] 1243.74195

25. Santoro, R; Elishakoff, I. Accuracy of the finite difference method in stochastic setting. J. Sound Vib.; 2006; 291, 1–2 pp. 275-284. [DOI: https://dx.doi.org/10.1016/j.jsv.2005.06.038] 1243.74195

26. Grigoriu, M., Khater, M.: Finite difference analysis of stochastic beams on elastic foundation, Report 85-7, Department of Structural Engineering, Cornell University, Ithaca (1985)

27. Seide, P. Accuracy of some numerical methods for column buckling. J .Eng. Mech. Div.; 1975; 101, 5 pp. 549-560. [DOI: https://dx.doi.org/10.1061/JMCEA3.0002050] 0312.73057

28. Thomée, V.: From finite differences to finite elements: a short history of numerical analysis of partial differential equations. In: Numerical Analysis: Historical Developments in the 20th Century, pp. 361–414, Amsterdam: Elsevier (2001)