Analysis of stochastically parameterized prestressed beams and frames

The stochastic behaviour of materials and loading is of great importance for buckling and analysis of prestressed beam and frame structures. In this paper, the stochastic stiffness and stress stiffness matrices are developed for stochastic analysis and buckling analysis. The bending rigidity and the buckling load are modelled as stochastic fields. The spectral decomposition known as the Karhunen–Loéve expansion has been used to expand the random fields. Using the Karhunen–Loéve expansion, the stiffness and stress stiffness closed-form matrices are formulated in terms of discrete parameters. The matrices are developed considering both classical Euler– Bernoulli theory and Timoshenko beam theory. Two case studies involving a pinned-pinned column and a frame structure is used to demonstrate the effectiveness of the proposed methods.


Introduction
In columns and frame structures prestressing can be of special interest.The influence of membrane forces on the lateral deflections can be considered by including the stress stiffening effect in the finite element (FE) analysis.Including the stress stiffening in the FE model also makes it possible to evaluate the buckling failure mechanism [1,2].Buckling failure in columns is characterized by large deflections perpendicular to the column axis when exposed to axial compression [2,3].Long slender beams loaded with axial compression are often more at risk than shorter thicker beams, as the critical buckling load will be reached before the critical yield load.
The critical Euler load, at which linear buckling occurs, was first described by Euler in 1744 [1].It is a well-known failure mode and methods for determining the buckling load have been developed intensively since the first publication [4,5].The linear buckling problem is an eigenvalue problem, which is relatively easy to solve in the case of simple structures [6][7][8].
It is often required to assess the uncertainties in engineering problems during the design phase.Uncertainties in loading, material properties and geometry or even human error can be included in the design by predetermined safety factors [2].These factors can be found in standards and guidelines such as the Eurocode [9] depending on the engineering problem at hand.Another method is to use probabilistic design, where the overall safety level (or probability of failure) is calculated [10].
The stochastic finite element method (SFEM) is one possible method to analyse the uncertainties in structural responses when spatial randomness is present in the structure [11][12][13].SFEM is an extension to the classic deterministic FEM, making it suitable also for problems that are otherwise difficult to assess analytically.SFEMs are well developed and have been applied to different types of practical engineering problems as described in the review paper by Arregui-Mena et al. [14].Different approaches can be used in the context of SFEM, including Monte Carlo simulation, perturbation approaches and Spectral SFEM.These methods have all been developed extensively over the years [11,12,15,16].
A great deal of research effort has been carried out in the field of stochastic finite element analysis (SFEA) of buckling.Lin and Kam [17] used the perturbation approach to formulate a stochastic FE (SFE) model for buckling of frames with uncertainty in geometric imperfections, material and sectional properties.The geometric imperfections were modelled using a Fourier series where the amplitudes were treated as random variables.Köylüoğlu et al. [18] predicted the mean and coefficient of variation of the stochastic buckling load using the SFEM for columns.The bending rigidity was modelled as a Gaussian field using three random variables by the weighted integral method [19], and Monte Carlo simulation was used for estimating the response variability.The critical buckling load was considered deterministic.Ramu and Ganesan [20] considered randomly fluctuating Young's modulus and stochastic axial loads in buckling of a column.The perturbation approach was used to express the eigenvalues and eigenvectors to determine the mean and covariances.Ramu and Ganesan [21] obtained the stochastic response of a column resting on a continuous Winkler foundation and discrete elastic supports.The Young's modulus, mass distribution and axial load were considered stochastic distributions.Vryzidis et al. [22], investigated the stochastic stability of steel tubes https://doi.org/10.1016/j.engstruct.2021.113312Received 15 January 2021; Received in revised form 10 September 2021; Accepted 29 September 2021 with random imperfections.The geometric imperfections in the steel tubes were modelled as a 2D univariate homogeneous stochastic field using a spectral representation method.A Monte Carlo based SFE formulation was implemented to investigate the response.A stochastic meshfree method for buckling analysis of columns was proposed by Gupta and Arun [23].The stochastic fields were discretized using the Karhunen-Loéve (KL) expansion and the method was used to estimate the response statistics of the buckling load considering different boundary conditions.Stochastic formulation of Timoshenko and Euler-Bernoulli theory using wavelet FE theory was described by Vadlamani and Arun [24].From this literature review, it can be concluded that significant effort has been put into stochastic analysis of beams and columns.However, to the best of the authors knowledge, no implementable closed-form and explicit solutions to both the stochastic stiffness matrix and the stochastic stress stiffness matrix considering random variability in the bending rigidity and prestress have been given previously.The main motivation of this paper is to develop explicit expressions to capture the effect of randomness in buckling and prestressed structures.Developing closed-form stiffness matrices makes it possible to easily include the stochastic stiffness in any general FE framework.Furthermore, closed-form stochastic matrices reduce the necessity of excessive Monte Carlo simulations compared to methods where the matrices are not directly provided.Adhikari and Friswell [25], developed closed-form solutions for the stochastic stiffness matrix considering classical Euler-Bernoulli beam theory.The KL expansion was used as basis for the random field expansion and up to 13 terms of the expansion were included.These matrices were used for model-updating of a stochastic undamped system.Adhikari and Friswell [25], did not consider the case of stochastic stress stiffening or stochastic buckling.Furthermore, only Euler-Bernoulli beam theory was used for the developed matrices excluding Timoshenko beam theory.
In the context of SFEA, the random fields can be discretized into random variables using the Karhunen-Loéve expansion [26,27].The KL expansion is often used for stochastic field representation in the context of SFEM [11,12,28].Generally, for complex random fields including more dimensions and certain correlation functions, numerical methods need to be used.But for 1D fields, analytical expressions for the KL expansion is given in the literature [11,25].
An additional aim of this paper is to develop a framework for predicting the stochastic response of beam buckling and prestressed frames, considering both the Euler-Bernoulli and Timoshenko beam theory.In the investigated literature, the stochastic analysis of beams and columns are often made only considering Euler-Bernoulli beam theory.Vadlamani and Arun [24], developed a stochastic formulation for Timoshenko beams using wavelet FE theory, but to the authors knowledge, no other effort has been made for including Timoshenko beam theory for SFEM of buckling and prestressed frames.Based on the literature review the following gaps in already published methods for stochastic buckling analysis have been identified: • Closed-form solutions to the stochastic stiffness and stress stiffness matrix for Euler-Bernoulli beams and Timoshenko beams have not been provided previously for use with classic SFEM.• No research effort has been carried out to use the KL expansion for field discretization for column buckling problems in the case of the regular FEM.The benefit of using the KL expansion is its simplicity and the presence of an analytical solution for 1D fields.• Previous works have mostly considered the classical Euler-Bernoulli beam theory.As influence from shear is relevant for shorter thicker beams, stochastic analysis of beams considering Timoshenko beam theory is required.
In this paper, the KL expansion is used to develop the closedform stochastic stiffness and stress stiffness matrix for beam and frame structures.The closed-form stochastic matrices are developed considering both classical Euler-Bernoulli beam theory and Timoshenko beam theory.The developed matrices are easy to implement in the regular FE formulations, making it possible to do SFEA without considerable effort.Furthermore, the closed-form matrices allows for efficient and fast pre-processing of the FE models, as the stochastic stiffness matrices are computationally fast to assemble.The developed framework is used on two numerical examples considering buckling of a beam structure and a prestressed frame structure.The outline of the paper is: In Section 2 the spectral decomposition of stochastic fields using the KL expansion is described and discussed for Gaussian fields.In Section 3 the KL expansion is used to expand the stiffness and stress stiffness matrices for an Euler-Bernoulli beam.This formulation is expanded to Timoshenko beams in Section 4. Based on the derived stochastic matrices buckling analysis of a beam structure and stochastic analysis of a prestressed frame structure is presented in Section 6.Finally, in Section 7, a set of conclusions is drawn.

Spectral decomposition of stochastic fields
In most static and dynamic problems some uncertainty, being in either geometry, loading or material properties can be expected.In these problems, the random parameters can be treated as random fields (stochastic fields) and SFEM can be used.Both dynamic and static problems considering randomness in stiffness, mass and damping, have previously been investigated inside the framework of SFEM [15,[28][29][30][31].In this paper, the stochastic response of stress stiffened beam and frame structures including critical Euler loads, are considered.Consider a random field (, ) with its covariance function   (  ,   ) defined in the space .  is a outcome from the random sample space , thus  ∈ .Each realization   corresponds to a single field (or function) (,   ).Generally, random fields are difficult to handle mathematically in the context of stochastic partial differential equations, such as the equations of motions [25].However, one possible method to handle this is by discretizing the random fields based on random variables.In this paper, the Karhunen-Loéve expansion is used for discretizing the random fields.The KL expansion given in Eq. ( 1), is a generalized Fourier type series based on the decomposition of the covariance function.
Where   () are uncorrelated random variables of some distribution type and  0 () is the deterministic part of the random field.In the rest of the paper (⋅) 0 refers to the deterministic part of (⋅).  and () are the eigenvalues and eigenfunctions of the integral equation given in Eq. (2).Eq. ( 2) is a Fredholm Integral of the second kind.
The KL series can be truncated to some finite number of terms, depending on the desired accuracy of the discretization.Generally, the Fredholm integral in Eq. ( 2) can be difficult to solve analytically, thus numerical integration methods are proposed [32,33].In this paper, the random fields considered are one dimensional.For one dimensional random fields, analytical solutions to the integral equation can be found in literature [11,25].The exponential autocovariance function is used in this paper and is expressed in Eq. (3).
Where  is the correlation length.The correlation length reflects the rate at which the correlation function decays between two points.If the correlation length is very large, the random field will behave as a single random variable system.If the correlation length is very small it results in a delta-correlated field.Exponentially decaying autocovariance is a popular choice when homogeneous random fields are used in stochastic analysis.A modified exponential autocovariance function [34], is sometimes preferred due to the non-differentiability of the exponential autocovariance function at its peak.In this paper, however, the correlation lengths are relatively long, reducing the error occurring due to the non-differentiability [11].
The random field  (, ) can be expanded as shown in Eq. ( 4), given that the mean is zero and that the domain is given in the interval Where analytical solutions have been derived for the eigenvalues and eigenfunctions [11,16,25].The eigenvalues and eigenfunctions for odd  are given as: where tan(  ) =    (5) and for even  they are given as: ,   () = sin(  ) Where  = 1∕.In practice the infinite series needs to be truncated to a finite number of terms.In [25], the number of terms to retain is described based on the 'amount of information' to keep.This is directly correlated to the eigenvalues,   .In Fig. 1, an example of the eigenfunctions scaled by the eigenvalues using Eqs.( 5) and ( 6) is plotted.As seen from Fig. 1, the amplitude of the eigenfunctions decrease with increasing index .Thus, the influence of the eigenfunctions in the discretized stochastic field decreases with increased number of KL terms.If one desire to keep 90% of the information in the KL expansion, the required number of KL terms,  can be determined as √   ∕ √  1 = 0.1 [25].The discretization is highly dependent on the correlation length .The greater the correlation length, the lesser number of KL terms is required as a high correlation length reduces the number of variables needed to describe the random field.In Fig. 2, the eigenvalues found using Eqs.( 5) and ( 6) are plotted for two different correlation lengths,  = ∕2 and  = ∕5, with  = 0.4.The number of KL terms required to obtain a 90% level of information is plotted.As seen from the plot a total number of 10 KL terms is required for  = ∕2, whereas 19 KL terms are required for  = ∕5.

Stochastic prestressed Euler-Bernoulli beams
The equation of motion for a static Euler-Bernoulli beam with random bending stiffness and axial load distribution is given in Eq. ( 7) [35].Where  () is the transverse displacement, (, ) is the flexural rigidity,  (, ) is axial load and () is the applied forcing.It is assumed that the bending stiffness, (, ) and axial force,  (, ) are random fields given in the form: where   are deterministic constants given as 0 <   ≪ 1 ( = 1, 2).
The stochastic fields   (, ) are taken to have zero mean, unit standard deviation and unit covariance.Therefore,   acts as strength parameters that quantify the amount of randomness in the random fields, that are used to model the distributed parameters of the system [25].When the random fields are Gaussian, the exponential decaying autocorrelation structure translates to both the random fields (, ) and  (, ).This is due to the fact that Gaussian random fields are invariant under linear translation.Randomness in the bending rigidity, (, ), can be caused by small variations in the geometry or Young's modulus.Axial load variation dominates certain engineering problems.For example, the self-weight of the system can be considered by a stochastic field.For reliability analysis considering the buckling phenomenon an additional stochastic  (, ) can also be included as an additional uncertainty in the model.The expression for (, ) in Eq. ( 8) can be used to formulate the stochastic element stiffness matrix as: Where () is the shape function for a Euler-Bernoulli beam.Following the notation from Adhikari and Friswell [25] the shape function can be defined as () =  (), where: Using the KL decomposition to expand the stochastic field in Eq. ( 10) we get: Where the deterministic part of Eq. ( 12) is given by: Extending Eq. ( 13) results in the deterministic stiffness matrix that can be found in most FE textbooks like [8]: The random part of Eq. ( 12) is given by Eq. ( 15): (15) Where   corresponds to the number of KL terms in the spectral decomposition.  are uncorrelated random Gaussian variables with zero mean and unit standard deviation.The constant matrices   are expressed by Eq. ( 16), for which Adhikari and Friswell [25] developed closed-form expressions.
Where   are the eigenfunctions from the solution to Fredholm Integral Eqs. ( 5) and ( 6).  is the starting position of the relevant element.
It should be noted that the KL expansion is used to discretize the stochastic field over the full length of the physical member (including all finite elements).Thus,   ensures that the correct part of the field is coupled to the correct finite element in the considered member.
Extending the work of Adhikari and Friswell [25], the stress stiffening matrix   , including the stochastic expression for  (, ) can be derived.The stress stiffening matrix is given as: Using the KL expansion on the stochastic field  2 (, ) we get: Where the deterministic part is given in Eq. ( 19) Which, when extended, results in the deterministic stress stiffening matrix: The random part of Eq. ( 18) is given as: Where   is the number of KL terms in the expansion of the random field of  (, ) and   () are uncorrelated Gaussian random variables with zero mean and unit standard deviation.The constant matrices of the stress stiffness matrix can be expressed as Eq. ( 22): Where the eigenfunctions from the KL expansion are defined in Eq. ( 5) and ( 6).Closed-form expressions of the matrices are derived in Appendix A.
After using the closed-form expressions to obtain the element stiffness   () and stress stiffness   () the usual methods can be used to assemble the global stiffness matrices.The global stiffness matrices will likewise be in the form of a deterministic part and a random part: where the deterministic parts,   and   are the general deterministic global stiffness matrix and global stress stiffness matrix, which can be found in FE textbooks like [8].The random parts can be assembled as: Where the element matrices   and   have been assembled into the global matrices   and    .The number of random variables   and   depends on the truncation of the KL expansion.

Stochastic prestressed Timoshenko beams
For shorter and thicker beams where the influence of shear deformation and rotational bending influences the deformation results, the Timoshenko beam theory [2,36,37] can be used.Similar to the stochastic matrices developed using the Euler-Bernoulli beam theory described above, the stochastic stiffness matrix and stochastic stress stiffness matrix, considering Timoshenko beam theory are derived.
The equation of motion for static Timoshenko beams with random bending rigidity and axial load distribution is given in Eq. ( 27) [35].

𝐸𝐼(𝑥, 𝜃) 𝑘𝐴𝐺 𝜕𝑝(𝑥) 𝜕𝑥
) Where  () is the transverse displacement, (, ) is the flexural rigidity,  (, ) is axial load and () is the applied forcing. is a constant known as the Timoshenko shear constant that depend on the cross section details.For rectangular sections  = 5∕6, can be used [38]. is the cross section area and  is the shear modulus.
It is assumed that the bending stiffness, (, ) and axial force,  (, ) are random fields given in the form: where   again are deterministic constants given as 0 <   ≪ 1 ( = 1, 2).
The stochastic fields   (, ) are taken to have zero mean, unit standard deviation and unit covariance.

Stochastic stiffness matrix
The expressions for (, ) can be used to formulate the stochastic element stiffness matrix: Where the contribution from bending   and shear   to the stiffness have been separated.The shear modulus can be found as  = ∕ (2(1 + )) for isotropic materials where  is Poisson's ratio.The shape functions follow the notation used for the Euler-Bernoulli beams with   () =     () and   () =     (): and Where  is defined as: gives the relative importance of the shear deformation to the bending deformation.Setting it to  = 0, results in the general Euler-Bernoulli shape functions and stiffness matrix. will be stochastic in nature as it consists of both (, ) and the shear modulus  which depends on .However, in the following derivations  is considered deterministic.Including the expression for the shear modulus for isotropic materials in Eq. ( 33) results in: As observed from Eq. ( 34), the resulting expression for  does not depend on the Young's modulus, .Only geometrical parameters  and  are required.Assuming a rectangular cross section with height ℎ and width , Eq. ( 34) can be rewritten into: As seen from Eq. ( 35),  is a function of ℎ.Thus, the variability of  is larger with increasing randomness in ℎ.However, the impact of  is small in the deterministic case, thus including randomness in  will lead to an even smaller contribution in the stiffness and stress stiffness matrix.
Using the KL decomposition to expand the stochastic field in Eq. ( 30), we get: Where the deterministic parts of Eq. ( 36) are given in Eqs.(37) and (38) and the random parts are given in Eqs.(39) and (40).(40) The deterministic parts of the equations can be given as the following deterministic stiffness matrix: In Eqs. ( 39) and ( 40), the constants   and   are the number of terms retained in the KL expansion and   and   are uncorrelated Gaussian random variables with zero mean and unit standard deviation.As both the bending contribution   () and the shear contribution   () in the stiffness matrix are derived from a random field of (, ), in all practical cases   =   and   =   .The constant matrices   and   can be expressed as: The eigenfunctions from the KL expansion are given in Eq. ( 5) and ( 6).Closed-form expressions are derived in Appendix B. Lastly the global stiffness matrix can be assembled using traditional methods in the form given in Eq. ( 23) to (26).

Stochastic stress stiffening matrix
The stress stiffening matrix, including the random field of the force  (, ) is given for a single element as: Where the shape function,   is the sum of the shape functions for bending and shear, expressed as   () =     (): Using the KL decomposition to expand the stochastic field in Eq. ( 44), we get: The deterministic and random part of the stress stiffening matrix are given as Eqs.( 47) and (48).
where the deterministic matrix can be expressed as: In Eq. ( 48),   is the number of terms retained in the KL expansion and   is the uncorrelated Gaussian random variables with zero mean and unit standard deviation.The constant matrix   can be expressed as: Where the eigenfunctions from the KL expansion are given in Eq. ( 5) and ( 6).Closed-form expressions are derived in Appendix B. Lastly the global stiffness matrix can be assembled using traditional methods, in the form given in Eq. ( 23) to (26).

Stochastic analysis
The expressions for the stochastic Euler-Bernoulli and Timoshenko beams are of special interest in two types of problems, namely classic beams and Euler buckling (eigenvalue buckling) analysis.In the context of SFEM, these two types of analysis are further described in this section.

Assumption of small deformations
The developed stochastic matrices in Sections 3 and 4 are based on the assumption of small deformations.This assumption is valid as long as the beams and columns are initially straight.However, if considerable initial bending is present or if the membrane deformation is highly coupled to the bending deformation, this assumption is no longer valid.In these cases, the buckling loads can be overestimated and non-linear analysis can be used instead [8].The developed matrices in this paper are, however, useful in the preliminary design stages for a given structure or for simple structures.Using the stochastic matrices, the stochastic output can be estimated and the sensitivity to bending rigidity or axial loading can be determined.

Stochastic response of prestressed beams
Stress stiffening is used to describe the lateral deflection of beams and columns influenced by membrane forces [8].When compressive membrane loads are added to a beam, the resistance to bending is decreased.The effect of membrane loading can be taken into account using the stress stiffening matrix, which is developed for stochastic beams in the above.The prestressing can be taken into account in the FE analysis by adding the stress stiffness matrix to the element stiffness matrix, as expressed in Eq. (51) [8].
The general FE approach for solving the displacements can be then be used: Where  is the externally applied force (deterministic) and  are the deflections.Stochastic analysis can be performed using classical Euler-Bernoulli beam theory or Timoshenko beam theory, using the developed closedform solutions provided in Appendices A and B. This results in the formal definition of the total stochastic stiffness given as: In the case where the influence of stress stiffening should be included in the model, the beam forces,  (, ), are required.The beam forces can be found using the general finite element approach and then be included in the assembly of the stress stiffness matrix.This means that a two-step calculation is necessary.In the first step, the beam forces can be found using the general FE formulation.From these forces, the stress stiffness matrix can be found and in the second step, the final FE analysis can be performed.Randomness in (, ) will not influence  (, ) in the case of simple beam buckling problems.External forces that are influenced by the beam properties, such as self-weight will, however, result in variations in  (, ).For the case of frame-structures, the distribution of forces will be altered based on the randomness in rigidity.For example, if the stiffness of one of the members changes, it will change the force distribution.The force in the frame itself will, however, remain constant.

Generalized random matrix eigenvalue problems
Buckling is the loss of stability in the equilibrium.For simple beam and columns (sometimes denoted beam-columns) the eigenbuckling, or bifurcation, has been described intensively for both analytical [1,3] and numerical solutions [6][7][8].Buckling can be considered a special case of stress stiffness analysis, where the membrane forces get so high that static stability is lost.
In the case of buckling, the governing equations of motion as given in Eq. ( 7) and ( 27) are considered.These equations neglect the inertia terms, making them static cases.In FE formulation, the corresponding deterministic buckling loads and buckling modes can be expressed as an eigenvalue problem in Eq. ( 54) [8].To solve the eigenvalue problem, the stress stiffness matrix,  , , is needed.The matrix can be determined with any arbitrary reference external load,   .The arbitrary load is used to calculate the forces  in the beams and columns, which are required to establish the stress stiffness matrix [8].

[ 𝑲 + 𝜆 𝑐𝑟,𝑗 𝑲 𝝈,𝒓𝒆𝒇
]   = 0 (54)  is the general static stiffness matrix.The matrices are provided for the stochastic case in the above.  are the eigenvalues (critical buckling multipliers) corresponding to the reference load and   are the eigenvectors (buckling modes).These should not be confused with the eigenvalues and eigenfunctions used with the KL expansion.Thus, for the traditional eigenbuckling problem, the resulting critical buckling load is predicted by the critical buckling factor   : For simple cases, where only a single axial load is applied to the structure, a unit load can be applied.This results in the   values corresponding directly to the critical load.Generally, only the first few eigenvalues and critical buckling loads are of interest.The higher modes are of less interest, as lower modes will result in buckling earlier than the higher modes.Using the KL expansions for the stress stiffness matrix and stiffness matrix, the random eigenvalue problem can be expressed as: Where each  , corresponds to a critical buckling multiplier for the reference load. , will be stochastic.
Eq. ( 56) can be solved using direct Monte Carlo analysis.Another possibility is to use the first-order perturbation approach for eigenvalue problems [39]: Where  0 and   are the deterministic eigenvalues and eigenvectors of the buckling problem.Eq. ( 57) is only valid, when the eigenvectors are scaled such that: Generally, the first-order perturbation approach is applicable for the first critical buckling loads.When higher-order buckling loads are sought, higher-order perturbation or other methods like the chaos expansion methods [40], should be used.
In the cases where the axial forces are not simple to estimate, a two-step calculation is required similar to the analysis described in Section 5.2.In the first step, the initial deflections in the structure are estimated and used to evaluate the element forces.The element forces are required to find the stress stiffness matrix.In the next step, the stochastic buckling equation Eq. ( 56) can be solved.When the element forces are found considering stochastic bending rigidity, (, ) the resulting forces will be stochastic.Thus, the stochastic stress stiffness matrix, which is built on the assumption that  (, ) can be directly found from the stochastic element forces.Additional randomness from self-weight or residual stresses can be included additionally by adding it to  (, ).

Numerical case studies
In this section, two numerical case studies are presented to estimate the response statistics using the developed stochastic matrices.In these case studies, the statistical moments are predicted using direct Monte Carlo analysis.In the case of buckling analysis, the perturbation approach has additionally been used for comparison.The matrices developed using the Timoshenko and classical Euler-Bernoulli theories are used and the output statistics have been compared.
In the first case study, the deflections in a prestressed redundant frame structure, that was investigated by [41], are evaluated.The analysis has been performed considering random variability in the bending rigidity.A deterministic preload has been added to the model.In the second case study, a simple pinned-pinned column is analysed by considering random variability in both the bending rigidity  and axial self-weight  , to show the methodology for buckling analysis.Results are obtained using both direct Monte Carlo sampling and the perturbation approach.

Case study of stochastic redundant frame structure
In this case study, a redundant frame structure [41], with variability in the bending rigidity () has been considered.The frame includes a bolting mechanism, which makes it possible to apply pretension in one of the diagonal frame members.The redundant structure is shown in Fig. 3, together with the corresponding FE model developed.The frame is built by rectangular cross-sections of 8 × 30 mm.In-plane bending is considered in the FE analysis.The two bottom edge nodes are fixed and simply supported respectively and two loads are applied as shown in Fig. 3(b).
It should be noted that two nodes are present in the frame on top of each other.This results in two diagonal frame members, which are not sharing a node at the centre, effectively uncoupling them.Two numerical cases have been investigated.In the first case, no pretension has been applied but stress stiffening effects have been included and in the second case, a pretension force of  = 500 N (compression) has been applied.In this numerical example, the pretension,  , is considered deterministic.For small loads, the frame with the shown tightening mechanism will be able to transfer both compression and tension loads.The unperturbed physical and geometric properties are given in Table 1.The bending rigidity is assumed modelled using a Gaussian distribution with a coefficient of variation of 8%.The correlation length is assumed to be  = 1∕5.Each of the frame members is modelled using their own KL expansion.This resembles the case where 6 individual members are used in the frame structure.Thus, zero correlation between the members is assumed.For the simulation of correlated random fields on multiple and non-convex domains, we refer to [42] for further details.The level of retained information has been set to 10%.This means that the longest frame members (the diagonals) require 21 KL terms, the horizontal frame members require 19 terms and the vertical members require 13 KL terms in the KL expansion.
The FE model is discretized using 30 elements for each frame members, resulting in a total element count of 180.This has been verified by a deterministic convergence study of the vertical deflection in the centre-node of the prestressed frame member, (highlighted in Fig. 3(b)).The convergence study is made on the non-prestressed model but considering stress stiffening effects.The mesh convergence study is shown in Fig. 4(a).If no stress stiffening effect was included, that is by excluding the stress stiffening matrix in Eq. ( 41), the deflection would converge faster.As the stochastic field is discretized through the FE mesh, five perturbations of (, ) are plotted in Fig. 4(b), for the left vertical frame member.As seen from the plot, the perturbations are rather fine with no steep changes, indicating a sufficiently refined mesh.The KL expansion is independent of the FE mesh, however, when applying the stochastic field to the FE model, the field has to be discretized to the mesh.Thus, a too course FE mesh will result in an unrealistic discretization of the stochastic field and thereby result in errors in the stochastic analysis.
The formulation given in Eq. ( 52) has been considered in a Monte Carlo calculation scheme where the total stochastic stiffness is found using Eq. ( 53).Both Timoshenko and Euler-Bernoulli beam theories have been applied.The analysis is considered in two steps.In the first step, the deflections and resulting axial forces in each element are found, based on the stochastic bending rigidity.The axial forces in each frame element are found using: Where  and  is the horizontal and lateral node deflections and  and  are the direction cosines in the frame element axis.In these calculations,  is considered deterministic but could be accounted for by the KL expansion given in the above sections.
It should be noted that the structure in this example is build up by frame elements as compared to the formulation given in Sections 3 and 4 for beam elements.This means that an additional degree of freedom (axial) has to be added and that rotation of the deterministic and stochastic stiffness matrices are required.This rotation can made using the traditional rotation matrix  , such that the rotated stiffness and stress stiffness matrix can be found by   =  ⊺  .The rotation matrix can be found in most books on FE analysis such as [8].The frame elements [8], have three DOF's, that is two translational DOF's and one rotational DOF, at each node.The axial degree of freedom is added to the stiffness matrices developed in Sections 3 and 4. As the axial degree of freedom is included in the stiffness matrix only and as it only depends on  and  it has been considered deterministic.
A total of 10,000 MC samples have been generated and used for the stochastic analysis.This has been verified by a convergence study, that showed that the mean and standard deviation of the output statistics converges around 6,000 to 10,0000 MC samples.To compare the results of the numerical study, the vertical deflection in the centre-node has been investigated.In Table 2, the results are given for the case where no pretension is added in the FE model and in Table 3, the results are given for the case where a pretension of  = 500 N is added.
It can be observed from the results from the stochastic analysis in Tables 2 and 3 that the Euler-Bernoulli and Timoshenko beam theories give nearly the same standard deviations and mean deflections.This is also observed in the cases with and without preloading.As the frame members are relatively thin and long, the difference between Euler-Bernoulli and Timoshenko beam theories is also expected to be small.The preload  = 500 N in the diagonal frame member results in smaller differences between deterministic and stochastic results, as seen in Table 3.As the preload is considered deterministic, it is expected that the deterministic and stochastic values will come closer.The beam forces in the diagonal preloaded frame member are approximately 65 kN, thus the additional preload of 500 N corresponds to 0.77% of the beam forces.This increase of force in the frame member effectively reduces the variability observed in output response.

Case study of stochastic buckling of a pinned-pinned column
In this case study, a simple pinned-pinned column with variability in the bending rigidity (, ) and the beam load  (, ) is considered for eigenbuckling.The variation in  (, ) is included in the case study for completeness.The column is shown in Fig. 5(a) and the unperturbed physical and geometric properties are given in Table 4.In the FE model, the mean of the force  (, ) is equal to 1.This results in the output critical load amplifier,   from the calculations, being equal to the stochastic load required for the column to buckle,   .The slenderness ratio, calculated as (  2 ) ∕ () is for the mean values calculated to 320.51.The slenderness ratio determines the effect of the Timoshenko beam formulation compared to the Euler-Bernoulli beam formulation.The lower the slenderness ratio, the higher the influence of shear  deformations on the results.From the formulation, it is seen that lower and thicker beams and columns result in lower ratios.From classical Euler-Bernoulli beam theory, the critical load for a pinned-pinned column like the one in Fig. 5(a), is 14.21 MN.
The column has been analysed using the Euler-Bernoulli and Timoshenko beam theories for comparison.The FE model of the column is developed using 40 elements.The convergence study of the critical buckling load using deterministic Timoshenko beam theory is shown in Fig. 5(b) and shows that the FE model converges by using 18 elements.However, the discretization of the stochastic field through the FE elements are rather coarse with 20 elements as seen from Fig. 6(a).Instead, 40 elements have been used, resulting in a finer discretization of the stochastic fields as seen in Fig. 6(b).
In the buckling calculations, it is assumed that both variational parameters can be modelled by a Gaussian random field.For the bending rigidity a 8% variation is assumed and a correlation length of  = ∕5.The critical buckling load is assumed to have a correlation length of  = ∕2 and a 5% variation.The KL expansion is truncated to obtain a 10% level of retained information.This corresponds to 19 KL terms for the rigidity  and 10 KL terms for the load  .This also corresponds to the eigenvalue plots shown in Fig. 2. As the load  (, ) only needs 10 KL terms, the KL discretization of the stiffness , has been used for the study of required elements shown in Fig. 6.Smaller correlation lengths results in a higher number of KL terms and in practice means that more fluctuations of the field can be expected.
To obtain the output statistics, both direct Monte Carlo simulations of the stochastic fields and the perturbation approach are used.The uncorrelated random Gaussian variables,   () for each KL expansions, are generated using a random number generator and have been used as input in the stochastic analysis.A total of 10,000 MC simulations for each random field have been generated.As multiple uncorrelated random variables are required in the KL distribution, the number of MC simulations is corrected for 19 and 10 KL terms respectively.Thus, for the bending rigidity where 19 KL terms are required a total of 190,000 random variables are generated.The eigenbuckling load multiplier is calculated from the eigenvalue problem given in Eq. ( 56), for each of the MC simulations.The critical load for obtaining buckling can then be determined as   =    .Only the first critical load corresponding to the first mode of flexure is estimated.From the critical loads, the statistics are estimated and provided in Table 5.The statistical moments using the perturbation approach is predicted using additional MC simulations and Eq.(57).Other methods for determining the statistical moments using the perturbation approach can be found in [43].It can be observed that the Euler-Bernoulli beam theory predicts the mean of the buckling load higher than the Timoshenko beam theory using both the Monte Carlo simulation method and the perturbation approach.This is also to be expected, as the Timoshenko beam theory will reduce the buckling load.It can also be seen that the two theories in combination with the perturbation method and the Monte Carlo simulation, results in almost the same mean and standard deviations.The highest difference between method is observed to be 2.25% for the mean value and 3.51% for the standard deviation.In Fig. 7, the probability density functions (PDF) are presented for the Euler-Bernoulli and Timoshenko     beam theories.As observed from Fig. 7 that the results from the Monte Carlo simulation and perturbation approach are very similar.It should be noted, that both the Monte Carlo approach and the perturbation approach and the beam theories rely on the KL expansion, which is truncated to a finite number of terms.Thus, both methods will result in approximated output statistics.
The output statistics of the buckling problem is sensitive to the boundary conditions.For example, changing the column in Fig. 5(a) to being fixed-free (cantilever beam), but keeping the rest of the assumptions, produces the results shown in Table 6.As observed from the results, the standard deviations considerably decrease and the relative differences between methods and theories likewise decrease.This indicates that the results are highly influenced by the boundary conditions.The PDFs in Fig. 8 are showing the same trend as in Fig. 7, as the results from the Monte Carlo simulation and perturbation approach are very close.So it can be concluded that both the Monte Carlo simulation and the perturbation approach give the same predictions.
In the above examples, the force  (, ) has been considered stochastic for completeness.A real life example could be in the case where self-weight is included in the FE model.Thus, considering the pinned-pinned beam in Fig. 5(a) the force through the beam could be modelled by: Where the density (, ) is considered a stochastic variable and  and  is the cross sectional area and gravitational constant, respectively. 1 is a deterministic force added on top of the beam.The result from the buckling analysis, the critical buckling factor,   , should be scaled with the reference loads, see Eq. ( 55) to obtain the critical loads.
In the following case study the deterministic buckling forces for the first buckling mode including self-weight has been applied as external loading  1 .The bending rigidity and density is considered stochastic following the properties given in Table 7.
Five perturbations of the stochastic force  (, ) is shown in Fig. 9.The force will increase stochastically from the top to the bottom of the beam, due to the stochastic density.
The resulting stochastic load factors,   from the analysis including self-weight is given in Table 8.As seen from the table, the results show similar tendencies for both theories and methods.It can be concluded that including the self-weight of the system results in stochastic critical load factors.

Conclusions
In this paper, closed-form solutions for the stochastic stiffness and stress stiffness matrix for beams and frames have been developed considering classical Euler-Bernoulli and Timoshenko beam theories.The developed matrices are based on the Karhunen-Loéve expansion and allow for simple implementation of additional KL terms, which results in more accurate stochastic analysis.The main motivation for developing the closed-form matrices arises from the fact that many engineering structures are at risk of local buckling.The proposed matrices include random variability in bending stiffness and prestress/buckling loads.The matrices can be applied to other random properties, following the methodology outlined in the paper.Furthermore, the assembly of the closed-form matrices is easy and corresponds to assembling general deterministic matrices in the FE formulation.The KL expansion is also independent of the number of elements, meaning that increased number of elements do not require additional random variables in the stochastic formulation.Thus, it can be easily implemented for random fields in stochastic buckling analysis of beams and frames by reducing the necessity for crude Monte Carlo sampling of each parameter for each element.The novel aspects in this paper include: a. Development of explicit closed-form solutions to the stochastic stiffness matrix and stochastic stress stiffness matrix considering variability in the bending rigidity and preload using the Karhunen-Loevé expansion.b.Development of the stochastic stiffness and stress stiffness matrices for Timoshenko beam theory.
where  = ∕2 and  corresponds to the full length of the beam.K  is expressed by: The elements of the symmetric matrix K  is then expressed as: For notational purposes the shape function   () is expressed as the product of   () =     () as shown in Eq. (31).  () is expressed as the product of   () =     () as shown in Eq. (32).The constant matrix given in Eq. ( 42) for the bending contribution is expressed as:

B.1. Contribution from bending -odd
Inserting the solution to the eigenfunctions from the KL expansion given in Eq. ( 5), the constant matrix for odd  can be expressed as: where  = ∕2 and  corresponds to the full length of the beam.K  can be expressed as: The elements in the symmetric matrix K  can then be expressed as:

B.2. Contribution from bending -even
For even  the bending contribution to the stiffness matrix can be expressed based on Eq. (6) The elements in the symmetric matrix K  can then be expressed as:

B.3. Contribution from shear -odd
The constant matrix given in Eq. ( 43) for the odd shear contribution to the stochastic stiffness can be expressed as: Inserting the solutions to the eigenfunctions from the KL expansion given in Eq. ( 5), the constant matrix for odd  is: The elements in the symmetric matrix K  can then be expressed as: The elements in the symmetric matrix K  can then be expressed as:

Appendix C. Closed-form expressions of stochastic stress stiffness matrix of Timoshenko beam
This appendix derives the closed-form expressions for the random part of the element stress stiffness matrix for Timoshenko beams.The matrix is defined as given in Eq. (44) as:

C.1. Odd
Inserting the solution to the eigenfunctions from the KL expansion given in Eq. ( 5), the constant matrix of odd  can be expressed as:

Fig. 1 .
Fig. 1.Eigenfunctions scaled with eigenvalues for the first 10 KL terms.An exponential covariance function and unit domains have been used.The correlation length is  = ∕5.

Fig. 2 .
Fig. 2. Number of KL terms required to obtain a level of retained information of 90% for different correlation lengths.

Fig. 6 .
Fig. 6.Example of perturbations of  considering various number of elements.

Fig. 8 .
Fig. 8. Probability density functions of critical buckling loads for fixed-free column.

Table 1
Material and geometric properties of the redundant frame structure.
Fig. 3. Frame structure and corresponding FE model.

Table 2
Vertical deflections at centre-node from stochastic analysis for the case without pretension added.

Table 3
Vertical deflections at centre-node from stochastic analysis for the case where pretension of  = 500 N is added.

Table 4
Material and geometric properties of pinned-pinned column.

Table 5
Critical buckling loads   from stochastic buckling analysis of pinned-pinned column.

Table 6
Critical buckling loads   from stochastic buckling analysis of cantilever beam.

Table 7
Material and geometric properties of pinned-pinned column exposed to self-weight.
a Critical buckling force using Euler-Bernoulli beam theory.b Critical buckling force using Timoshenko beam theory.

Table 8
Critical buckling factors   from stochastic buckling analysis of pinned-pinned column.