Probabilistic Engineering Mechanics The exact element stiffness matrices of stochastically parametered beams

Stiffness matrices of beams with stochastic distributed parameters modelled by random fields are considered. In stochastic finite element analysis, deterministic shape functions are traditionally employed to derive stiffness matrices using the variational principle. Such matrices are not exact because the deterministic shape functions are not derived from the exact solution of the governing stochastic partial differential equation with the relevant boundary conditions. This paper proposes an analytical method based on Castigliano’s approach for a beam element with general spatially varying parameters. This gives the exact and a simple closed-form expression of the stiffness matrix in terms of certain integrals of the spatially varying function. The expressions are valid for any integrable random fields. It is shown that the exact element stiffness matrix of a stochastically parametered beam can be expressed by three basic random variables. Analytical expressions of the random variables and their associated coefficient matrices are derived for two cases: when the bending rigidity is a random field and when the bending flexibility is a random field. It is theoretically proved that the conventional stochastic element stiffness matrix is a first-order perturbation approximation to the exact expression. A sampling method to obtain the basic random variables using the Karhunen–Loève expansion is proposed. Results from the exact stiffness matrices are compared with the approximate conventional stiffness matrix. Gaussian and uniform random fields with different correlation lengths are used to illustrate the numerical results. The exact closed-form analytical expression of the element stiffness matrix derived here can be used for benchmarking future numerical methods.


Introduction
In the context of continuum modelling of complex engineering systems, probabilistic uncertainties can be modelled as random variables or random fields [1]. When random field models are used to represent spatially varying uncertain quantities, the Stochastic Finite Element Method (SFEM) [2][3][4] is the systematic approach available for such problems. The implementation of the stochastic finite element method has two major steps, namely, (a) the incorporation of the random fields within the scope of the finite element modelling to derive the stochastic governing equations and the boundary conditions, and (b) to solve the stochastic equations using uncertainty propagation approaches. The second step has received significant attention as this is often the computationally expensive part. Several reduced-order computationally efficient methods, such as first-and second-order perturbation methods [5,6], Neumann expansion method [7,8], polynomial chaos method [2] and spectral function method [9], have been developed. The first step, although not always the computationally most expensive step, governs the accuracy of the overall results. This paper is about this first step larger elements, this fact may not be true. In this paper, this aspect of stochastic modelling is investigated using Euler-Bernoulli beams as an example.
Mechanics of beams with stochastic field parameters is a classical topic. In one of the earliest works, Vanmarcke and Grigoriu [17] discussed stochastic finite element analysis of beams. Deodatis [18] and Deodatis and Shinozuka [19] considered static finite element analysis of Euler-Bernoulli beams with random field properties. They have also proposed a formulation based on stochastic shape functions. Chang and Chang [20] considered stochastic dynamic finite element analysis of a beam on a random foundation. Elishakoff et al. [21,22] proposed exact solutions for bending problems of a beam with spatially stochastic properties. Ren et al. [23] considered variational principle based finite element method for stochastic beams. The book by Elishakoff and Ren [24] present a comprehensive account of finite element for stochastic mechanics problems when the underlying variability is large. For dynamic problems, several authors have proposed spectral formulations in the frequency domain for beams with random field properties [25][26][27]. In [28], a doubly spectral stochastic finite element method was developed for Euler-Bernoulli beams based on the dynamic shape functions of the underlying deterministic problem. More recently, Larsen et al. [29] developed a static shape function based stochastic finite element method for Euler-Bernoulli and Timoshenko beams for buckling analysis. From this brief review of literature, it is clear that an exact stiffness matrix of a beam with parameters modelled as random fields is not available in literature in closed-form. Direct numerical simulation is possible to obtain the exact stiffness matrix within a Monte Carlo simulation framework. However, such approaches lose the elegance and simplicity of closed-form stiffness matrix expressions available for deterministic beams.
Motivated by such needs, this paper aims to develop Euler-Bernoulli beams' exact element stiffness matrix with spatially randomly varying parameters. Unlike the conventional variational formulation used in the stochastic finite element formulation, a direct force-displacement approach based on Castigliano's method has been developed. As a consequence, the use of the finite element shape functions is not necessary. The board aims of the proposed formulation are: • to develop the exact and simple closed-form expression of the element stiffness matrix of Euler-Bernoulli beams with general spatially varying parameters • to consider random fields describing both the bending rigidity and flexibility in a unified manner • to have the ability to include Gaussian and non-Gaussian random fields within the same analytical framework • to keep the number of random variables at the minimum to reduce further numerical works (such as the structural reliability analysis) It will be shown that only three random variables are necessary to include the effect of random fields exactly in the stiffness matrix, irrespective of the nature of the random fields. The outline of the paper is as follows. Overview of the conventional stochastic stiffness matrix for Euler-Bernoulli beams derived using the cubic shape functions is given in Section 2. In Section 3, Castigliano's approach for a beam element with general spatially varying parameters is developed. This gives the exact closed-form expression of the stiffness matrix in terms of certain integrals of the spatially varying function. In Section 4, the special case when the bending flexibility is a random field is considered. It is proved that the exact stiffness matrix can be expressed by three basic random variables, which are linear functions of the underlying random field. The case, case when the bending rigidity is a random field is discussed in Section 5. It is proved that the exact stiffness matrix can be expressed by three basic random variables, which are nonlinear functions of the underlying random field. In Section 6, the theoretical developments are numerically illustrated. The statistical properties of the basic random variables, elements of the stiffness matrix and the response of a cantilever are illustrated. In the Appendix, a sampling method to obtain the basic random variables using the Karhunen-Loève expansion is proposed. This is achieved through an explicit closed-form expression of the covariance matrix derived by evaluating the underlying integrals analytically. Finally, some key conclusions are drawn in Section 7.

Overview of the conventional stochastic stiffness matrix
In this section we briefly review the conventional approach to obtain the stiffness matrix of a beam with random field properties. The equation governing the deformation of an Euler-Bernoulli beam [30] of length with random bending rigidity can be expressed as In this equation ( ) is the transverse bending displacement, ( , ) is the bending rigidity, is the length along the beam and ( ) is the applied distributed forcing. Here denotes an element of the (random) sample space so that ∈ . It is assumed that the bending rigidity is a random field of the form In Fig. 1, a beam with a random bending rigidity is shown. In an alternative way, the flexibility function can also be modelled as a random field. In that case one has The subscript 0 indicates the mean values, 0 < ≪ 1 ( = 1,2) are deterministic constants (also known as the strength parameters) and the random fields ( , ) are taken to have zero mean, unit standard deviation and covariance ( 1 − 2 ) (with a stationary random field assumption). Since, ( , ) and 1∕ ( , ) are strictly positive functions ∀ , , ( , ) ( = 1,2) are required to satisfy the conditions This requirement, strictly speaking, rules out the use of Gaussian models for these random fields. However, for small , it is expected that Gaussian models still can be used if the primary interest of the analysis is to estimate the first few response moments and not the response behaviour near tails of the probability distributions. The analytical formulation to be developed in this paper is not dependent on the Gaussian or stationary assumptions of the random fields. However, these assumptions often simplify subsequent numerical calculations and interpretation of the results. The beam element in Fig. 1 has two nodes and four degrees of freedom when only the bending deformation is considered. The displacement field within the element is expressed by cubic shape functions [30,31] for the classical finite element analysis and they are given by where the constant matrix and the vector function are Note that these cubic shape functions arise from the exact solution of the governing differential equation of the deterministic system underpinning Eq. (1) with relevant boundary conditions. Employing these Fig. 1. Depiction of an Euler-Bernoulli beam with bending rigidity modelled as a random field, ( , ). The beam element has four degrees of freedom, corresponding to transverse deformation and rotation at the two ends.
shape functions in conjunction with the conventional variational formulation [30,31], the stiffness matrix of the general beam element can be obtained aŝ In the above (•) ′′ stands for double derivative with respect to . Separating the integral, one obtainŝ where the deterministic part is given bŷ The integrand matrix in the above equation is given by Evaluating the integral in Eq. (9) and simplifying, the deterministic part of the stiffness matrix is obtained aŝ This is the conventional stiffness matrix of an Euler-Bernoulli beam [30][31][32]. The random part in Eq. (8) is given bŷ The random variables ( ) are defined as These random variables appear as power integrals of the random field (equivalent to moments over the spatial domain). They were termed as 'weighted integral' in [33,34] and was introduced for stochastic static finite element analysis of Euler-Bernoulli beams by Deodatis [18] and Deodatis and Shinozuka [19]. Later Manohar and Adhikari [35] extended this approach for dynamic problems when the weighted integral appear in terms of integrals over transcendental functions. The random variables ( ) are correlated random variables as they are derived from the same random field. In this paper, the random variables ( ), and other similar random variables to be introduced later, are called basic random variables. These are a fundamental minimum number of random variables necessary to represent the total randomness of the system modelled using random fields. The random fields appearing in the definition of the weighted integral in Eq. (22) are general as long as the physical condition in Eq. (4) is satisfied. In the special case, if the random field is assumed to be stationary (homogeneous), continuous and Gaussian, then the random variables ( ) in Eq. (13) will also be Gaussian. For such cases, the random field can be discretised into random variables using the Karhunen-Loève (KL) expansion [36][37][38][39][40]. The discretised variables can, in turn, be used to obtain statistical moments and consequently probability density functions of the random valuable ( ). For a certain form of the autocorrelation function, such as the exponential function, the Karhunen-Loève expansion can be obtained in closed-form [2,41]. Using such expansions, the random valuable ( ) are explicitly obtained in Appendix.
The expression of the random part of the stiffness matrix can be rewritten such that it is a linear combination of the basic random variables ( ). From Eq. (12) one haŝ The constant matriceŝ, = 1, 2, 3 are obtained aŝ Similar expressions were obtained by Deodatis [18]. Some of the key observations of practical interests from this formulation are: • The random and the deterministic parts of the stiffness matrix are linearly separated. • The random part of the stiffness matrix is a linear superposition of three correlated random variables only. This is irrespective of the nature of the underlying random field (homogeneous/inhomogeneous or Gaussian/non-Gaussian). • If the random field is Gaussian, the three random variables appearing in the expression of the stiffness matrix will also be Gaussian. Fig. 2. Schematic diagram of a general two-noded beam member with a variable cross-section. The quantities 1 and 1 are the shear force and bending moment, respectively at node 1. The quantities 2 , and 2 are the shear force and bending moment, respectively at node 2.
Although the random part of the stiffness matrix is obtained analytically in closed-form, it is not exact. This is because the shape functions used to derive them is not obtained from the exact solution of the stochastic ordinary differential equation Eq. (1). An alternative approach is proposed to derive the exact stiffness matrix in the next section.

The Castigliano's approach for a beam element with general spatial variabilities
This section proposes the procedure to obtain the closed-form expression of a beam's stiffness matrix with variable cross-section. Castigliano's method is utilised to derive the exact stiffness matrix in [42,43]. To explain the essential equations of Castigliano's approach, the generalised force and displacements are shown in Fig. 2. To obtain the stiffness matrix, a flexibility-based approach is considered. The matrix form of the force-displacement relation for a general beam can be expressed as = (17) Here and are the generalised force and displacement vectors of the following form and and is the 4 × 4 stiffness matrix. The quantities and ( = 1, 2) are the nodal bending displacements and rotations of the beam crosssection, respectively. To obtain the stiffness coefficients, the flexibility approach is considered followed by the equilibrium conditions. To derive the force-displacement relationship for node 1, node 2 is kept fixed (all two degrees of freedom are restrained), and generalised forces 1 and 1 are applied on node 1.
Only bending deflections are considered in the formulation. So, the strain energy is expressed as The internal moment is given by = 1 − 1 . The flexibility relationship of the beam is given by The coefficients s are defined as The stiffness relationship can be found by inverting the flexibility equation (21) as where, 1 = 1 3 − 2 2 . In a similar way, the direct flexibility matrix for point 2 can be obtained. For that, we have to fix point 1 and apply Castigliano's 2nd theorem after putting the internal moment equation = 2 ( − ) − 2 in Eq. (20). Therefore, one obtains Here s are given by the following integrals The stiffness relationship for node 2 can be found by inverting the flexibility relation Eq. (24) as where, 2 = 1 3 − 2 2 . The coefficients and are as follows and The determinant Now, considering Eqs. (23), (26) and (27) The stiffness relationship for the whole beam is given below The stiffness terms corresponding to the coupling between nodes 1 and 2 can be found considering the moment and force equilibrium as For the force equilibrium one obtains 1 = − 2 . Considering Eqs. (37) and (38) in force equilibrium equation one can obtain the relationship between some of the coupling terms. Those are given below as ( 11 + 31 ) 1 + ( 12 + 32 ) 1 + ( 13 + 33 ) 2 + ( 14 + 34 ) 2 = 0 (41) Equating the coefficients of the left hand and right hand side we get and The above relationships can be written as Now, the stiffness matrix become Considering the moment expressions of moments and shear the moment equilibrium equation 1 + 2 + 2 = 0 becomes 1 + 2 + 2 = 21 Equating the coefficients of the left hand and right hand side we obtain The complete bending stiffness matrix is therefore obtained by combining the above relations as This is the most general, concise and exact expression of the bending stiffness matrix with variable bending flexibility. The only restriction is that 1∕ ( ) is finitely integrable as given by Eq. (22). In the special case when ( ) is constant with respect to , the integrals s become Substituting these in Eq. (50), it can be seen that this stiffness matrix reduces exactly to the conventional expression in Eq. (11). Therefore, Eq. (50) can be viewed as the generalisation of the classical beam stiffness matrix with constant cross-section to a variable cross-section. For frames and other built-up structures, this element matrix can be employed with transformation matrices and assembled in the usual procedure used in the finite element method. Next, this general expression is applied when the flexibility and rigidity of the beam are modelled using random fields.

Bending flexibility is a random field
The expression of the exact stiffness matrix with arbitrary variable cross-section is given by Eq. (50). There are three integrals, namely For the random field 2 ( , ), we define the random variables ( ) as Using these definitions, the can be related to the random variables ( ) through the following relationships Substituting these, in the expression of the exact stiffness matrix in (50), after some algebraic simplifications, one obtains The functional dependence with is used to show the quantifies which are random. The expressions of 1 ( ) and the two matrices in Eq. (55) are given by and The constant matrices , = 1, 2, 3 are derived as From Eq. (55) it can be seen that when the beam is deterministic, that is, when 2 = 0, it reduces to the classical deterministic expression given by (11). For the random case, unlike the approximate analysis in Section 2, the stiffness matrix, in general, cannot be separated into a deterministic and random part due to the nonlinear nature of 1 ( ) in the denominator. Some of the key observations from the exact analysis when the bending flexibility is a random field are: 1. The random and the deterministic parts of the stiffness matrix cannot be separated in a linear manner. 2. In the limit, when the bending flexibility is deterministic, the random stiffness matrix reduces to the classical deterministic stiffness matrix. 3. Unlike the conventional analysis in Section 2, the random variables do not appear as a linear superposition. The expression of the stiffness matrix is a nonlinear function of the three random variables. 4. Similar to the conventional analysis in Section 2, the random variables are linear functions of the underlying random field. 5. If the random field is Gaussian, the three random variables appearing in the expression of the stiffness matrix will also be Gaussian.
In the next section, we consider the case when bending rigidity is random.

Exact expressions
When the bending rigidity is a random field, ( , ) is given by Eq. (2). It is necessary to evaluate the integrals 1 , 2 and 3 given in Eqs. (27)-(29) as they appear in the expression of the exact stiffness matrix in Eq. (50). Using the expression of the random bending rigidity function in Eq. (2), the integrals can be expressed as By rearranging the denominator one has Using this, the integrals in Eq. (60) can be rewritten as In the above, the new random field ( , ) is defined as a nonlinear function of the random field 1 ( , ) through For the random field ( , ), we define the random variables ( ) as Using these definitions, the integrals s can be related to the random variables ( ) through the following relationships similar to what discussed in the previous section Substituting these, in the expression of the exact stiffness matrix in (50), one obtains Although the expressions of the stiffness matrix are similar to the case of random flexibility discussed in the previous section, there is a crucial difference. The random field ( , ) in Eq. (63) is a nonlinear function of the random field 1 ( , ). As a result, the random variables ( ) will be non-Gaussian even when the random field 1 ( , ) is Gaussian. Therefore, the stiffness matrix in Eq. (66) is in general a non-linear function of non-Gaussian random variables.

Approximate analysis
In many real-life engineering problems, one often observes small variabilities in the material and geometric properties due to stricter quality control. As the random fields have zero mean and unit standard deviation, the randomness is quantified by the strength parameters , = 1, 2. It is assumed that Using Taylor series expansion, we can establish From the inequality of the norms (Cauchy-Schwarz inequality) we have This is because, due to the unit standard deviation, | 1 ( , )| = 1 when a probabilistic average is taken. Therefore, the higher order terms above 2 1 in Eq. (68) can be neglected. Considering the definition of the random field 2 ( , ) in Eq. (3), one can deduce that when 1 = 2 For notational convenience we denote As the expressions derived in Section 4 are relatively simple, we employ them to derive the first-order perturbation expressions. From Eq. (56), the inverse of 1 ( ) can be approximated as Substituting this in Eq. (55) and keeping only the first-order terms in one can express the stiffness matrix as Here the deterministic part 0 is the same as in (11) and the random part can be expressed as a linear combination of the random variables ( ) in (13) as Comparing the expression of the random part of the stiffness matrix in Eq. (75) with the conventional expression in Eq. (14), it is observed that they are identical. This is a remarkable coincidence. It proves that the conventional expression of the random element stiffness matrix is equivalent to a first-order perturbation approximation of the exact expression. This exciting fact has not been established yet as the exact general expressions were not available in the literature. Based on this mathematical derivation, we have this fundamental remark in stochastic finite element analysis: Remark 1. The conventional element stiffness matrix of an Euler-Bernoulli beam with a random field bending rigidity derived using the cubic shape functions is a first-order perturbation approximation of the exact element stiffness matrix.
An important question that can be asked at this point is whether this result is general. That is, can we draw the same conclusion, for example, for Timoshenko beams, plates, shells and other 3D elements. If such a generality can be established, then it will have a huge impact in quantifying errors in the stochastic finite element analysis as the conventional stiffness matrix approach discussed in Section 2 is usually used in practice.

Numerical results and discussions
Theoretical derivations and underlying concepts proposed in the previous sections will be numerically demonstrated here. A beam of length = 1.0 m and a nominal cross-section of 0.03 × 0.003 m 2 are considered. The material property assumed is that of aluminium, that is, = 69 × 10 9 N/m 2 . Therefore, the nominal bending rigidity of the beam is 0 = 4.66 N m 2 . The autocorrelation function of the random fields to be employed in this study is the exponential function as given in Eq. (A.85). For numerical calculations, we discretise the random fields using the Karhunen-Loève (KL) expansion. Further discussions on KL expansion can be found, for example, in Refs. [2,12,40,41,[44][45][46][47]. Karhunen-Loève expansion of a random field is an infinite series. We truncate this series as in Eq. (A.89) with a finite value of . The number of terms in the KL expansion depends on the correlation length and amount of information to be retained in the series. In general, the smaller the correlation length, the more terms are necessary for the KL expansion.
To obtain the three different basic random variables ( ), ( ) and ( ), for = 1, 2, 3, given by Eqs. (13), (53) and (64), it is not necessary to use the KL expansion. We use KL expansion only for computational purposes. One can use direct simulations methods [48] to generate the random fields and then numerically integrate to obtain the random variables.
Once the truncated Karhunen-Loève (KL) series (A.89) is established, it can be employed for different basis random variables , = 1, … , . The following two cases are considered for numerical illustrations: • Case 1: All the basis random variables are independent and identically distributed (i.i.d.) Gaussian random variables with zero mean and unit standard deviation, that is ∼ (0, 1), ∀ . The correlation length in the exponential autocorrelation function is considered to be ∕10. For this case, = 56 number of random variables are used. This captures 94% of the information in the KL expansion.
• Case 2: All the basis random variables are i.i.d. uniform random variables with zero mean and unit standard deviation, that is ∼ (0, 1), ∀ . A larger correlation length of ∕2 is considered. For this case, = 18 number of random variables are used, capturing 95% of the information in the KL expansion.
These two cases are designed to give a broader perspective on the nature of the numerical results. The studies in the rest of this section consider these two cases.

Case 1: Gaussian random variables
The three different basic random variables ( ), ( ) and ( ), for = 1, 2, 3, given by Eqs. (13), (53) and (64), underpin the expressions of the three stiffness matrices developed in the paper. Random parts of the stiffness matrices, given by Eqs. (14), (58) and (66) utilise these three basic random variables in unique ways. It is therefore of primary interest to understand the nature of the basic random variables. The random variables ( ) and ( ) are statistically identical in nature provided the same random field is used for the bending rigidity and flexibility. Therefore, we only consider the random variables ( ) and ( ), which are derived using different approaches. In Fig. 3, the probability density functions of the random variables ( ) and ( ), for = 1, 2, 3 are shown. The pdfs are obtained with 10,000 samples with Gaussian random variables (case 1). The random variables ( ) are sampled using the analytical covariance matrix given in Eq. (A.102). The samples of ( ) are calculated by performing the integration in Eq. (64) numerically. Note that ( ) change with . However, negligible difference was observed in our numerical calculations for different values of . Therefore, results for = 0.10 are shown as representative values in Fig. 3.
From the three probability density function plots in Fig. 3, it can be observed that ( ) are ( ) are close for all = 1, 2, 3. Their nature is close to Gaussian random variables as all the basis random variables are Gaussian. The covariance matrices of ( ) are ( ) are calculated as The covariance matrix of ( ) is obtained directly from the closed-form expressions in Eq. (A.102). The covariance matrix of ( ) is calculated using the numerical integration in Eq. (64) with 10,000 Monte Carlo samples. The numerical integral for each sample involved 56 functions in the numerator and denominator. The resulting covariance matrix as shown in Eq. (78) is close that of ( ).

Case 2: Uniform random variables
In a strict sense, Karhunen-Loève expansion is applicable to Gaussian random fields only. This is because the sum of non-Gaussian random variables is, in general, not Gaussian. However, due to the central limit theorem (see, for example, [39]), if many terms are used in the KL expansion with a non-Gaussian random variable, the resulting random field tends to be a Gaussian random field. The advantage of using uniform random variables in KL expansion is that, unlike the previous case, parameters can be selected such that it always produces positive random fields. This is a necessary requirement for random fields to be used to model bending rigidity, which must always remain positive for stable structures. In Fig. 4, the probability density functions of the random variables ( ) and ( ), for = 1, 2, 3 are shown. The pdfs are obtained with 10,000 samples with 18 uniform random variables (case 2). The random variables ( ) are sampled using the analytical expression Eq. (A.92). Unlike the Gaussian case, they were not sampled from the analytical covariance matrix given in Eq. (A.102). However, numerical results show that they are very close. The samples of ( ) are obtained by performing the integration in Eq. (64) numerically with = 0.10. From the three probability density function plots in Fig. 4, we observe that ( ) are ( ) are close for all = 1, 2, 3. Their nature is, however, significantly different from Gaussian random variables as all the basis random variables are uniform. This is a major difference from  The covariance matrix of ( ) is calculated using the numerical integration in Eq. (64) with 10,000 Monte Carlo samples. The numerical integral for each sample involved 18 functions in the numerator and denominator. The resulting covariance matrix as shown in Eq. (79) is close that of ( ). We only show the second-order properties of the random variables. As seen in Fig. 4, the random variables are highly non-Gaussian in nature. Therefore, higher-order moments may be useful to characterise them.

Characteristics of the stiffness matrix coefficients
Once the basic random variables are obtained, they contribute to the random stiffness matrix. Three approaches to obtain the stiffness matrix have been proposed. They are (a) conventional approach with random bending rigidity given by Eq. (8), (b) the exact approach with random bending flexibility given by Eq. (55), and (c) the exact approach with random bending rigidity given by Eq. (66). Three different notations, namelŷ( ), ( ) and̃( ), have been used to distinguish the three expressions. The three random element stiffness matrices use the random variables ( ), ( ) and ( ), for = 1, 2, 3, given by Eqs. (13), (53) and (64), respectively. In this section, the statistical nature of the element of the stiffness matrices obtained using these three different approaches will be discussed.
For easier comparative analysis, the random variables are divided by their corresponding deterministic values given by Eq. (11).

Case 1: Gaussian random variables
In Fig. 5, the four elements of the stiffness matrix obtained using the three different approaches are shown. Probability density functions calculated using a 10,000-sample Monte Carlo simulation are compared. Twenty percent randomness, that is, = 0.2 is considered. Recall that the conventional approach, the exact random bending flexibility approach and the exact random bending rigidity approach use the basic random variables ( ), ( ) and ( ), for = 1, 2, 3 respectively. They, in turn, are obtained using 56 basis i.i.d. Gaussian random variables as discussed in the previous section. Due to the use of Gaussian random variables, from Fig. 5 it is observed that the four random stiffness coefficients are Gaussian in nature. They have a nominal mean of 1 (as they are normalised by the corresponding deterministic values) and similar standard deviations. Scrutinising further, it can be seen that the results from the two exact approaches are closer to each other than the conventional method. This is expected as it was mathematically proved that the conventional method is a first-order perturbation approximation to the exact approach.

Case 2: Uniform random variables
We consider the same four elements of the stiffness matrix as before, namely 11 ( ), 12 ( ), 22 ( ) and 24 ( ). The probability density functions of these selected normalised stiffness coefficients are shown in Fig. 6 for = 0.2. These pdfs are obtained using a 10,000-sample Monte Carlo simulation with uniform random variables resulting the basic random variables ( ), ( ) and ( ), for = 1, 2, 3. The nature of the pdfs is very different from the Gaussian case discussed in the previous subsection. The pdfs of all the four normalised coefficients are centred about one but have a flatter region, resembling a uniform distribution. The results from the two exact approaches are closer to Fig. 5. The probability density function (pdf) of four selected normalised stiffness coefficients. The pdfs are obtained with 10,000 samples with Gaussian random variables (case 1) and the value of = 0.2 has been used. The random stiffness coefficients are normalised with respect to their deterministic values given by Eq. (11). The conventional approach, the exact random bending flexibility approach and the exact random bending rigidity approach use the basic random variables ( ), ( ) and ( ), for = 1, 2, 3 respectively. each other and significantly different from the conventional method. This implies that for non-Gaussian basis random variables, relatively more error is shown by the approximate conventional method.

Response analysis of a cantilever beam
In Fig. 7, a cantilever beam with a random bending rigidity is shown. A force is applied at the tip which results a deflection under the force. As the bending rigidity is a random field, the tip deflection will be a random variable. The tip deflection of the nominal beam is given by For the random case, eliminating the degrees of freedom at the fixed end, we have The random tip deflection ( ) is obtained by solving this equation. Direct Monte Carlo Simulation is used to obtain the samples of ( ) are the results for the two parametric cases are described next.

Case 1: Gaussian random variables
The three different basic random variables ( ), ( ) and ( ), for = 1, 2, 3, given by Eqs. (13), (53) and (64) are used to obtain three different stiffness matrices as explained before. For this case, 56 i.i.d. Gaussian random variables are simulated to generate the results. In Table 1, we show mean and standard deviation of the normalised tip deflection using 10,000 samples. The normalised tip deflection of the cantilever is obtained as Results for four values of the strength parameter is shown in Table 1.
The value of standard deviation increases with increasing values of as expected. The difference between the exact and approximate results with random rigidity increases with higher values of . Note that uncertainty propagated exactly in all three methods using the direct Monte Carlo simulation. Therefore, the differences observed arise purely due to the formulation methods and not due to the uncertainty propagation methods.
In Fig. 8, the probability density functions of the normalised deflection are shown for the four values of . The nature of the pdfs is close to Gaussian distributions. The results from the two exact methods are closer to each other and different from the approximate conventional method.

Case 2: Uniform random variables
The three different basic random variables ( ), ( ) and ( ), for = 1, 2, 3, are obtained with 18 i.i.d. uniform random variables. The nature of these random variables and the resulting stiffness coefficients was discussed before. A key observation was that the probability density functions of the basic random variables and the stiffness coefficients were non-Gaussian in nature. Here we examine if the same trend is followed to the response also.
In Table 2 we show mean and standard deviation of the normalised tip deflection. Results for four values of are shown, and they are calculated using 10,000 samples as before. The difference between the exact and the approximate method is more compared to the Gaussian case. In Fig. 9, the probability density functions of the normalised deflection are shown for the four values of . The results for the two exact methods are closer to each other. However, significant differences between approximate (conventional) and exact results are observed, especially for larger values of . This shows the impact of the formulation of the stochastic problem. There can be errors in the response statistics even if an exact simulation-based uncertainty propagation method is employed.

Conclusions
This paper addressed a fundamental issue in the stochastic finite element analysis, that is, the incorporation of random fields within the scope of the finite element modelling. A four-node Euler-Bernoulli beam element was used as an example to illustrate three distinct possibilities for stochastic finite element modelling. The three approaches to obtain the element stiffness matrix discussed are: • Conventional with stochastic rigidity: The random and the deterministic parts of the stiffness matrix become linearly separated. The random part of the stiffness matrix is a linear superposition of three random variables, which, in turn, are linear functions of the underlying random field. • Exact with stochastic flexibility: The random and the deterministic parts of the stiffness matrix are not linearly separable. The random part of the stiffness matrix is a nonlinear function of three random variables, which, in turn, are linear functions of the underlying random field. • Exact with stochastic rigidity: The random and the deterministic parts of the stiffness matrix are not linearly separable. The random part of the stiffness matrix is a nonlinear function of three random variables, which, in turn, are also nonlinear functions of the underlying random field.
The last two approaches are possible due to the derivation of the exact stiffness matrix using Castigliano's approach for a beam element with a general spatially varying parameter. The resulting closed-form expression is simple and valid for any finitely integrable function representing the spatial bending flexibility. In the special case, when the bending flexibility is a constant, the general stiffness matrix reduces exactly to the conventional expression. Each of the derived stiffness matrices is expressed as a combination of three basic random variables and constant matrices. Closed-form expressions of the constant matrices have been derived for all three cases. The basic random variables appear as power integrals of the underlying random fields. In general, they should be obtained numerically. However, if the random field is a stationary and Gaussian random field, then Karhunen-Loève expansion can be used to represent them. In that case, the integrals leading to the basic random variables have been evaluated in closed form. This gives an analytical expression of the covariance matrix, from which the samples can be generated easily. In the numerical analysis, the statistical properties of the basic random variables, elements of the stiffness matrix and the response of a cantilever are compared between the three proposed approaches. Noticeable differences were observed when uniform random variables were used in the simulation of the underlying random fields.
One of the important theoretical results that emerged from this study is that the first-order perturbation approximation of the exact element stiffness matrix with stochastic rigidity is the same as the conventional stiffness matrix with stochastic rigidity derived using the usual cubic shape functions. It proved that the conventional expression of the random element stiffness matrix is a first-order perturbation approximation of the exact expression. This result puts a clear perspective on the accuracy of conventional stochastic finite element analysis. Future research is needed to verify whether this result is general and applicable to other problems such as stochastic Timoshenko beams, plates, shells and general 3D elements. Only static problem is considered in this work. Future work is necessary to consider dynamic problems and the derivation of the exact, consistent mass matrix for stochastic problems.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix. Closed-form expressions of the basic random variables using Karhunen-Loève expansion
Three different basic random variables, involving power integrals of the underlying random fields, are introduced for the bending deformation analysis of beams. They are ( ), ( ) and ( ), for = 1, 2, 3 given by Eqs. (13), (53) and (64), The random variables ( ) and ( ) involve linear functions of the underlying random field. For this case, closed-form expressions can be obtained analytically as explained in this section.

A.1. Karhunen-Loève expansion
Suppose ( , ) is a random field with a covariance function ( 1 , 2 ) defined in a space . Here denotes an element of the (random) sample space so that ∈ . Mathematically and numerically, it is challenging to deal with random fields directly in the equations of motion, which are often expressed by partial differential equations. For this reason, it is required to discretise a random field in terms of random variables. Once this is done, then a wide range of mathematical and numerical techniques can be used to solve the resulting discrete stochastic differential equations. Among the many discretisation techniques, the spectral decomposition of random fields using the Karhunen-Loève expansion turns out to be very useful in practice. In this paper, this approach has been applied to model updating.
Since the covariance function is finite, symmetric and positive definite, it can be represented by a spectral decomposition. Using this spectral decomposition, the random field ( , ) can be expressed in a generalised Fourier type of series as The spectral decomposition in Eq. (A.83), which discretises a random field into random variables, is known as the Karhunen-Loève expansion.
The series in Eq. (A.83) can be ordered in a decreasing series so that it can be truncated after a finite number of terms with the desired accuracy. Fukunaga [38] and Papoulis and Pillai [39], and the references therein give further discussions on the Karhunen-Loève expansion.
In this paper one dimensional systems are considered. To demonstrate the approach a Gaussian random field with an exponentially decaying autocorrelation function is considered. Such a model is representative of many physical systems and closed form expressions for the Karhunen-Loève expansion may be obtained. The autocorrelation function between points 1 and 2 can be expressed as Here the constant is known as the correlation length, and it plays an important role in the description of a random field. If the correlation length is very small, then the random field becomes close to a delta-correlated field, often known as white noise. The random field effectively becomes a random variable if the correlation length is very large compared to the domain under consideration. Assuming the mean is zero, then the underlying random field ( , ) can be expanded using the Karhunen-Loève expansion [38,39] in the interval − ≤ ≤ as