Uncertain Dynamic Characteristic Analysis for Structures with Spatially Dependent Random System Parameters

This work presents a robust non-deterministic free vibration analysis for engineering structures with random field parameters in the frame of stochastic finite element method. For this, considering the randomness and spatial correlation of structural physical parameters, a parameter setting model based on random field theory is proposed to represent the random uncertainty of parameters, and the stochastic dynamic characteristics of different structural systems are then analyzed by incorporating the presented parameter setting model with finite element method. First, Gauss random field theory is used to describe the uncertainty of structural material parameters, the random parameters are then characterized as the standard deviation and correlation length of the random field, and the random field parameters are then discretized with the Karhunen–Loeve expansion method. Moreover, based on the discretized random parameters and finite element method, structural dynamic characteristics analysis is addressed, and the probability distribution density function of the random natural frequency is estimated based on multi-dimensional kernel density estimation method. Finally, the random field parameters of the structures are quantified by using the maximum likelihood estimation method to verify the effectiveness of the proposed method and the applicability of the constructed model. The results indicate that (1) for the perspective of maximum likelihood estimation, the parameter setting at the maximum value point is highly similar to the input parameters; (2) the random field considering more parameters reflects a more realistic structure.


Introduction
Structural dynamic characteristics, including the natural frequency and natural mode, as a crucial indicator for the vibrational properties of engineering structures [1], have been widely realized and studied for many years. With the help of finite element method (FEM), the structural dynamic characteristic can be adequately dealt with by investigating the generalized eigenvalue and eigenvector problems [2][3][4]. For example, Gorman and Yu [4] reviewed the method of superposition in vibration analysis of plates and shells, especially focusing on the Gorman method for accurate establishment of eigenvalues and mode shapes in free vibration analysis of rectangular plates. Although the system variables of the concerned structures are broadly accounted for as deterministic, it has been demonstrated that the fluctuation, i.e., uncertainties, of these parameters inevitably and inherently correlated to the structural modelling and analysis process [5][6][7]. The complexity of the actual structural material properties and various random errors during the manufacturing process will result in uncertainty of the structural parameters, such as vibration of the machine tool, random variation in the temperature during processing, etc., which will cause uncertainty among a group of structural components with the same nominal size that are manufactured with the same material and the same processing method [8] and ultimately lead to random fluctuation of material properties around the mean value and a certain correlation between the fluctuation and machining dimension direction [9,10]. In addition, a group of structural members with the same nominal size, because of the uncertainty of their material parameters, may have similar but different dynamic characteristics [11]. The existence of such uncertainties intrinsically has an influence on the believability of the analyzing results of the structural dynamic behaviors [12][13][14]. Hence, it is urgent to develop an uncertain free vibration analysis framework for more effective and meaningful estimation on the structural dynamic characteristics.
In general, uncertain dynamic characteristic analysis is implemented with probabilistic/stochastic approaches, which are based on the theory of probability or statistics [15][16][17][18][19]. Wan [15] used low-order statistical moments to adopt to characterize the uncertainty of modal frequencies of two bridges with assumed normally and uniformly distributed parameters. Liu et al [16] presented a probabilistic boundary element method for analysis of the statistics of structural eigenvalues and eigenvectors with random shape parameters. Consequently, in most of the literature listed above, the relevant uncertainties of the structural parameters are modelled as random discrete variables with predefined statistical information, such as mean values and variance. Over the past decade, stochastic finite element methods (SFEM) have been paid much attention and have been applied to structural analysis of static responses [20,21] or dynamic responses analysis [22][23][24], in which structural response problems were addressed by SFEM by incorporating probabilistic strategies within the FEM. Numerous computational procedures have been developed for solving the random static problems [25,26], as well as other engineering applications, including reliability problems [27]. Normally, there are two categories to implement SFEM: simulation approaches (e.g., the Monte-Carlo method), which are capable of offering the probabilistic features of the concerned structural responses based on the statistics of samples obtained from the simulation [28][29][30]; and non-simulation methods, which approach the statistical characteristics of the structural outputs by carrying out various numerical methods [31]. Development of SFEM in structural engineering, however, has not yet deeply and adequately extended to eigenvalue problems despite their importance in many applications, including the dynamic response of structures.
Moreover, admitting the universal application of SFEM, the creditability of such stochastic processes is conditional to the availability of the statistical information of the concerned uncertainties in practical engineering applications [32]. Especially, the research on characterization of structural uncertainty and modification of structural parameters mainly focuses on expansion of some structural or material parameters, such as elastic modulus, the moment of inertia, thickness, etc., and inputting the mean values and mean variances of these discrete random variables into the structural system for uncertainty analysis. Such constructed models are not very complex and comparatively easy to resolve in most cases but cannot accurately reflect the real uncertainty existing in input parameters and output outcomes of the actual structural system, although various effective methods can be used for stochastic finite element analysis, as listed above. In fact, there are many factors that affect the uncertainty of structural material parameters, resulting in the random distribution of material parameters in the structural space. For example, manufacturing processes can easily lead to spatial variations in the load and material properties, such as moduli and density. When rolling steel plates, the runout of the rolling head presents a trigonometric function law, which will inevitably cause uneven thickness of the structure and then variation in random physical parameters with structural spatial sizes. With the robust progress of uncertainty analysis, SFEM has been escalated with consideration of the spatial dependency of uncertain system parameters by incorporating the theory of random field with structural static analysis [6,10,33,34] or dynamic analysis of simple one-dimensional random field [35,36]. Moreover, there is still a lack of verification methods for SFEM with random field. In most cases, the Monte-Carlo method is used for such a purpose, but it definitely and inevitably has the disadvantage of absolute dependence on the amount of samples [37][38][39][40]. Hence, it is very crucial and urgent to build an effective analysis model of free vibration for the structure with random field and to develop a new verification method regarding the presented model.
Aiming at the uncertainty of structural material parameters, a method used of structural dynamic characteristics analysis and the corresponding verification are presented in the context of stochastic finite element method by incorporating the theory of random field with the finite element method. First, the uncertainty of structural material parameters is represented by the parameter setting model of random field theory, and the uncertainty of structural parameters is quantified with the random field model based on Gauss kernel function. Then, the simulation and discretization of the parameter setting model of random field are implemented with Karhunen-Loeve expansion method, and the structural dynamic characteristics are analyzed in the frame of finite element method, followed by acquisition of the probability distribution density function of the natural frequency by using the multi-dimensional kernel density estimation method. Afterwards, the input parameters of the model are quantified and verified by the maximum likelihood estimation method after comparing the experimental results with the simulation results. Finally, two examples are, respectively, used as one-dimensional and two-dimensional cases of random fields to validate the applicability and effectiveness of the proposed method.

Representation of the Uncertainty with Random Field
The errors of machining, heat treatment, and material itself may cause uncertainty of the structural system. These errors are usually small and independent. According to central limit theorem, the distribution of many independent and small random variables follow a Gaussian distribution. Actually, the assumption of Gaussian distribution is easy to calculate and the corresponding problems can be solved. In this work, Gauss random field model is used to describe the uncertainty of material parameters of structural system.

Gauss Random Field Model
Gauss random field has two characteristics: (1) its mathematical expectation µ and variance σ 2 are constants independent of position coordinates; i.e., m ω (x i ) = µ, D ω (x i ) = σ 2 , whereby ω(x i ) is a random number and x i represents a point in space; (2) its autocovariance function is uniquely related to the relative position distance of two points in the random field but not to the absolute position coordinate of two points; i.e., autocovariance C ω (x i , x i + τ) = C ω (τ) = σ 2 ρ(τ), whereby τ is the relative distance between two points, and ρ(τ) is the autocorrelation function of the random field. In addition, another important parameter of Gaussian random field is the correlation distance L, which indicates that the parameters within the correlation distance have obvious correlation. The key to establish a random field is to construct its covariance matrix.
In the framework of finite element, the continuous Gauss random field needs to be discretized into random variable vector for the following structural dynamic analysis. There are several commonly used discretization methods for Gauss random field, i.e., spectrum representation, Karhunen-Loeve expansion (K-L expansion), and so on.

Karhunen-Loeve Expansion
Karhunen-Loeve expansion has been widely applied to the continuous process [31]. Essentially, a random field is decomposed into a series of uncorrelated random variables and certain coefficients, such as eigenfunctions and eigenvalues, by using the K-L expansion. K-L expansion has the following advantages: it has the characteristics of mean square convergence for any type of random field; compared with other discretization methods, when the finite terms of expansion are the same, K-L expansion has the minimum mean square error [6]. In the K-L expansion, a random field H(x, θ) can be expanded into a group of countable and orthogonal random variables; i.e., H(x, θ) can be expanded into a combination of the random scalars ξ n (θ). The specific development process is as follows [1]: whereby x = (x, y) is the coordinate of a point in space, E(x) is the mathematical expectation, θ is a random event, ξ n (θ)(n = 0, 1, . . . , m) form a Gaussian random sequence with zero mean and they are uncorrelated to each other, and λ n and f n (x) are eigenvalues and eigenfunctions of the autocovariance matrix of random field C(x 1 , x 2 ) .
The key of K-L expansion for a random field is to obtain the eigenvalue and eigenfunction of autocovariance matrix C(x 1 , x 2 ) . Because the covariance matrix of a random field is defined in the regular geometric space domain, the eigenvalues and eigenvectors can be easily obtained. The detailed solution process of autocorrelation matrix C(x 1 , x 2 ) is as follows [5].
whereby the autocovariance function C(x 1 , x 2 ) is bounded, symmetric, and positive definite, which ensures that the eigenvalues and eigenfunctions have the following properties: (1) the set of eigenfunctions f i (x) is orthogonal and complete; (2) for each eigenvalue λ k , there are at most a limited number of linearly independent eigenfunctions; (3) there is at most one countable infinite set of eigenvalues; (4) all of the eigenvalues are positive real numbers; (5) the autocovariance function C(x 1 , x 2 ) can be decomposed into the following forms: In the case of one-dimensional (1D) random field, the following autocorrelation kernel function can be used for solving the eigenfunction f n (x) in Equation (2) C(x 1 , whereby the correlation length of random field L also reflects the attenuation degree of correlation between two points x 1 and x 2 . Therefore, C(x 1 , x 2 ) is a function of variable |x 1 − x 2 | and the parameter L; the integral area D in Equation (2) is a real number interval in the case of one-dimensional random field, and it can be taken as D = [−a, a]. Therefore, Equation (2) can be converted into: Equation (5) is further expanded to obtain: Finding the first derivative of Equation (6) to x 1 yields Finding the derivative of Equation (7) to x 1 once again and substituting x 1 with x, then the differential equation in general form can be obtained as follows Let ω 2 = − 2 L + λ n L 2 λ n ; Equation (8) can be converted into To solve Equation (9), its boundary conditions should be found at first. Substituting x 1 = −a and x 1 = a into Equations (6) and (7), respectively, and rearranging them, the following boundary conditions can be obtained In summary, Equation (2) has been transformed into a general differential Equation (9) with its boundary conditions Equations (10) and (11), and the general solution of Equation (9) is By substituting Equation (12) into boundary conditions and then rearranging the equations, this yields: Based on Equation (13), to obtain the solution of differential Equation (9), namely, to obtain Equation (12), the following conditions must be met: Solving Equation (14) and expressing the solution of the first equation as ω # n , ω # n should be within the interval n 2 − 1 2 π a , n 2 · π a , (n = 1, 3, . . . , m − 1), whereby m is an even number and it indicates the truncated number during the K-L expansion of Equation (1). Expressing the solution of the second equation of Equation (14) as ω * n , ω * n is within the interval n 2 − 1 2 π a , n 2 · π a , (n = 2, 4, . . . , m). The corresponding eigenfunctions are the following Equations (15) and (16) Furthermore, based on the transformation relation ω 2 = − 2 L + λ n L 2 /λ n , the corresponding eigenvalues to ω # n and ω * n can be derived as Based on Equations (14)- (16), take the Karhunen-Loeve expansion of beams for example. Set the beam length l = 1m, a = 0.5m and the correlation length L = 0.3m. When the number of truncations is taken as n = 1, 2, 3, 4, 5, 6, the eigenfunctions f n (x)  Figure 1. It can be seen directly that the eigenfunctions are actually a series of trigonometric functions, and their periods and amplitudes are related to the values of the chosen order n; that is, the larger the value of n, the smaller the period of f n (x). Based on Equations (17) and (18), when l = 5m, a = 0.5m, and L = 0.2m, 0.3m, 0.5m, 1m, 2m, 5m, the eigenvalues are shown in the right part of Figure 1. It can be seen from Figure 1 that the random field H(x, θ) is composed of eigenvalues λ n and eigenfunctions f n (x) based on Equation (1); the larger the value of L in the right part is, the larger the value of low-order eigenvalues, such as λ 2 , and the larger the proportion of low-order eigenfunctions, such as f 2 (x) in H(x, θ), as shown in the left part, so the fluctuation of random field is more gentle.   After the eigenfunctions and eigenvalues of 1D random field are obtained, the autocorrelation function can be obtained based on Equation (3); i.e., C(x 1 , When the length of 1D random field is 1m and the correlation length L is 0.3 m, C(x 1 , x 2 ) is displayed in Figure 2. Figure . When the length of 1D random field is 1m and the co relation length L is 0.3 m, ( ) , C x x is displayed in Figure 2. Figure     It can be seen from Figure 2 that ( ) 1 2 , C x x is more and more close to its autocorr lation kernel function with the increasing m ; that is, the more accurate the autocorrel It can be seen from Figure 2 that C(x 1 , x 2 ) is more and more close to its autocorrelation kernel function with the increasing m; that is, the more accurate the autocorrelation function is, the more precise the random field simulation is. Therefore, in the following examples in Section 4, the truncation number m is selected as 12 with a synthetical consideration of simulation accuracy and calculation workload.
For a two-dimensional (2D) Gaussian random field, its eigenvalues and eigenfunctions can be expressed by the product of the eigenvalues and eigenfunctions of two 1D random fields [3] as follows Substituting Equations (19) and (20) into Equation (1), the K-L expansion of 2D random field can be obtained

Structural Parameters Uncertainty Characterizing and Quantification
Based on the proposed Gaussian random field model, the uncertainty of structural parameters will be characterized and quantified in this section. In the framework of finite element, the Gauss random field of structural parameters is discretized into every grid element, and the random dynamic characteristics of the structure are then calculated. Nonparametric estimation of the structural dynamic characteristics will then be implemented by using the kernel density estimation method so as to obtain the distribution function characteristics of the structural random dynamic characteristics. Finally, the distribution characteristics of the output responses from the test and simulation will be compared with each other, and the distribution parameters of the model will be quantified and verification carried out with maximum likelihood estimation.

The Analysis of Structural Dynamic Characteristics with Random Field
Considering that the Young's modulus of the structure is a Gaussian random field, the random field is discretized based on Equation (1) and then substituted into the following element stiffness matrix and element consistent mass matrix (22) [ whereby [B] is the geometric matrix, [N] is the shape function matrix; Ω is the integral domain, Ω is the unit length for the bar element and beam element, and Ω is the unit area for the plate and shell elements; ρ is the material density; [D] is the elastic matrix shown in Equation (24) for the beam element, I z is the inertia moment of the beam; [D] is the bending stiffness matrix in Equation (25) for the thin plate element, µ is the Poisson's ratio of the material.
[D] = EI z (24) The total stiffness matrix [K] and the total mass matrix [M] of the structural system can be obtained by assembling element stiffness matrix and the element mass matrix and then substituting the boundary conditions of nodes into them. Generally, the vibration that causes damage to the structural system is low-frequency vibration, so, in the subsequent analysis, only the low-order natural frequency of the system is considered in this work, and the matrix iteration method is used for solution so as to quickly obtain the low-order natural frequency of the system and its corresponding vibration mode [6].
Suppose that the structural system represented by the stiffness matrix [K] and the mass matrix [M] is a positive system with n degrees of freedom and the free vibration equation of the structural system described with the flexibility matrix is as follows: ..
Let the solution of Equation (26) (27) can be converted into For the first-order eigenvalue λ 1 , the relation holds; i.e., [S]{X} 1 = λ 1 {X} 1 . Based on this relationship, the matrix iteration method is used to carry out iteration, and the maximal eigenvalue λ 1 and the corresponding eigenvector {X} 1 can be obtained. The specific iteration steps are as follows To solve the second-order and higher-order eigenvalues and eigenvectors, the dynamic matrix [S] needs cleaning; that is, [S] needs modifying by using the orthogonality of the main mode and the components of the first r-order main modes in [S] should be cleared so as to obtain the r + 1 order iterative dynamic matrix. According to the projection theorem of functional theory, the cleaning matrix can be obtained as: whereby Let the dynamic matrix before and after cleaning be [S] r and [S] r+1 ; the specific process of cleaning is as follows [10]: The r + 1 order eigenvalue of the system λ r+1 and its corresponding eigenvector {X} r+1 can be obtained by the operation of cleaning for [S] r+1 mentioned above and iterative calculation, and the r + 1 order natural frequency of the system can also be obtained, i.e., f r+1 = 1

Multidimensional Kernel Density Estimation
Maximum likelihood estimation provides a method to evaluate model parameters with given observation data; i.e., by observing the results of many tests, the parameter can be found that can make the probability of sample occurrence the maximum. Since the distribution characteristics of the output response of the model are unknown, it is necessary to make nonparametric estimation for the probability distribution of the output response. Herein, nonparametric estimation for the probability distribution of the output response is implemented with the multi-dimensional kernel density estimation method.
When the structural parameter is random field, the natural frequency f r (r = 1, 2, . . .) of the model is also random variable. Suppose that natural frequencies are independent and identically distributed random variables so that the multidimensional kernel density of its distribution density function is estimated as [8]: . . , f rm s j are the sample data of f r , f r is the mean value of sam- is the multivariate kernel function; h k is the window width matrix, the selecting principle of h k is to minimize the mean square error of calculation results; the kernel density estimation points q = (q 1 , q 2 , . . . , q k , . . . , q m ), q is a vector that is the independent variables of kernel density estimation function; herein, q is the vector . . , f rm s j .

Maximum Likelihood Estimation
The maximum likelihood estimation is implemented aiming at the output response of the structural model, and the parameter satisfying the probability distribution density is estimated. The parameter with the greatest possibility θ * is taken as the estimation value of the real parameter θ. Since the log likelihood function is easier to calculate, the kernel density estimation function in Equation (31) is then expressed aŝ whereby the parameter setting θ = (σ, L), σ is the mean variance of Gaussian random field, and L is the correlation length of random field; simulation output response Q mod is the output response of the model with parameter setting θ; , f exp is the output responses of tests, and num means the number of test points. According to the assumption of kernel function of multivariate Gaussian distribution [8], the window width is taken as The Gaussian kernel function is chosen as The maximum likelihood estimation method is used to quantify distribution parameters of test samples; i.e., assuming a group of parameters θ, the output responses of parameter samples are taken as the model samples. The probability distribution density function of the output responses is then obtained by using the kernel density estimation aiming at output response values based on Equation (31), and the maximum likelihood estimation functionL(θ) can be computed by inserting the output responses from both test samples and model samples into Equation (32), and, finally, parameter θ * can be estimated, which corresponds to the maximal value ofL(θ) so as to verify whether it is the input parameter of test samples θ. ( )

I-Beam with One-Dimensional Random Field
The maximum likelihood estimation method is used to quantify distribution eters of test samples; i.e., assuming a group of parameters θ , the output respo parameter samples are taken as the model samples. The probability distribution function of the output responses is then obtained by using the kernel density est aiming at output response values based on Equation (31), and the maximum lik estimation function  ( ) L θ can be computed by inserting the output responses fr test samples and model samples into Equation (32), and, finally, parameter * θ estimated, which corresponds to the maximal value of  ( ) L θ so as to verify whe the input parameter of test samplesθ .  Suppose that Young's modulus of material E is random field and its mean ( ) 210GPa E X  =  , and the covariance function is ( )

I-Beam with One-Dimensional Random Field
Then, pansion of E in every element from Equation (1) is (35), the random field E is represented with its mean vari and correlation length L; i.e., the parameter setting E θ is the output paramete samples and model samples, In the following, the random fiel Suppose that Young's modulus of material E is random field and its mean value is E(X) = 210 GPa, and the covariance function is C(x 1 , x 2 ) = e −|x 1 −x 2 |/L . Then, K-L expansion of E in every element from Equation (1) is From Equation (35), the random field E is represented with its mean variance σ and correlation length L; i.e., the parameter setting θ E is the output parameter of test samples and model samples, θ E = (σ E , L E ). In the following, the random field E is quantified with numerical tests. Based on Section 3.1, the first eight-order natural frequencies can be solved for every I-beam, and the natural frequency matrix can then be formed as Given an input parameter θ 0 = (σ 0 , L 0 ), and natural frequency matrix of 500 beams In the following, the kernel density estimation and maximum likelihood estimation will be implemented aiming at natural frequencies, and the identifiability of parameter setting θ = (σ, L) and the effectiveness of random field model will be verified.

Test Data Analysis
For four parameter settings of Young's modulus θ E1 = (1GPa, 100mm), θ E2 = (5GPa, 100mm), θ E3 = (1GPa, 1000mm), and θ E4 = (5GPa, 1000mm), the random process distributions of E on 10 beams are demonstrated in Figure 4, in which every curve denotes the value of E in the beam and reflects the randomness of values of E on every point of the beam. From Figure 4, it can be seen that values of the random field E fluctuate around its mean value 210 GPa; the fluctuation along the longitudinal axis reflects the magnitude of variance and the variation of variance with the axial dimension of the beam; just as the meaning of parameter setting σ E and L E , the influence of σ E on the random distribution of E is much greater than that of L E . In the case of the same parameter setting used, 10 curves in each figure, i.e., random distribution of E, are very different, but the general distribution law is similar in each figure. Moreover, it can be observed from Figure 4 that L E directly affects the randomness distribution of E on every element in the beam; that is, the larger the correlation length is, the more acute the value fluctuation of E along the direction of beam length.   Based on the conclusions of Figure 4, the parameter setting θ E = (5GPa, 500mm) will be used in the subsequent analysis.

One-Dimensional Kernel Density Estimation
In order to verify the applicability of the parameter setting θ E = (5GPa, 500mm), the test samples f First, take L E0 = 500mm and σ E ∈ {σ 1 , σ 2 , σ 3 , . . . , σ 10 } = {1, 2, 3, . . . , 10}GPa, which form 10 parameter settings θ E . For 100 beams corresponding to each θ E , f 2 is solved and the model sample { f 2 } mod 1×100 is formed. Then, take 100 beams with parameter settings σ E0 = 5GPa and L E0 = 500mm and f 2 of 100 beams construct test samples { f 2 } exp 1×100 . Figure 5a displays the one-dimensional kernel density estimation of f 2 calculated based on { f 2 } mod 1×100 and { f 2 } exp 1×100 . From Figure 5a, it can be seen that distributions of natural frequencies of model samples become more and more concentrated with the decreasing σ E , and the distributions are very close to each other when σ E of both test samples and model samples is 5 GPa.    Figure 5 shows that the influence of σ E on the fluctuation of kernel density is greater than that of L E .
In the same way, considering that the Young's modulus E and mass density ρ are random fields varying with the spatial size of beam, the parameter settings are taken as θ = (5, 500) and 100 beams are used. The 1D kernel density of f 2 is estimated again, where σ E = 5GPa and σ ρ = 5kg/m 3 . The left part is the 1D kernel density distribution function estimated of f 2 when L ρ = 500mm but L E changes, and the right part is the 1D kernel density estimated of f 2 when L E = 500mm but L ρ varies. When L ρ = L E = 500mm, in Figure 5a, σ E changes but σ ρ = 5kg/m 3 , and Figure 5b is the result estimated when σ E = 5GPa but σ ρ changes.
It can be seen from Figure 6 that the influence of σ E and σ ρ on the kernel density distribution function is significantly greater than that of L ρ and L E ; when other parameters are fixed but σ E and σ ρ change, respectively, the influence of σ E on the kernel density distribution function is greater than that of σ ρ ; when other parameters are fixed but L ρ and L E change, respectively, the influence of L E on the kernel density distribution function of f 2 is greater. In general, the random field E has a greater influence on the kernel density distribution function of f 2 than the random field ρ. eters are fixed but E σ and ρ σ change, respectively, the influence of E σ on the kernel density distribution function is greater than that of ρ σ ; when other parameters are fixed but L ρ and E L change, respectively, the influence of E L on the kernel density distribution function of 2 f is greater. In general, the random field E has a greater influence on the kernel density distribution function of 2 f than the random field ρ .

Multidimensional Kernel Density Estimation and Maximum Likelihood Estimation
In order to accurately verify the validity of the random field model and the identifiability of the parameter settings, it is necessary to further estimate the multi-dimensional kernel density probability distribution function of the structure output responses, and then to carry out the maximum likelihood estimation to obtain the parameter setting corresponding to the maximal value of likelihood function. Because the values of probability density function of some points are small, in order to avoid the calculation value of likelihood function being 0, it is necessary to take the logarithm of the value of probability density function first and then implement summation. The obtained log likelihood function is illustrated in Figure 7.

Multidimensional Kernel Density Estimation and Maximum Likelihood Estimation
In order to accurately verify the validity of the random field model and the identifiability of the parameter settings, it is necessary to further estimate the multi-dimensional kernel density probability distribution function of the structure output responses, and then to carry out the maximum likelihood estimation to obtain the parameter setting corresponding to the maximal value of likelihood function.
When considering random field E, based on 100 model matrix samples [f] mod 8×100 from 1 × 10 4 beams and test samples [f] exp 8×500 from 500 beams, the 8-dimensional kernel density of the first 8-order natural frequencies of the beam is estimated, and then the parameter θ E = (σ E , L E ) is estimated by the maximum likelihood based on Equation (32), in which Because the values of probability density function of some points are small, in order to avoid the calculation value of likelihood function being 0, it is necessary to take the logarithm of the value of probability density function first and then implement summation. The obtained log likelihood function is illustrated in Figure 7.   Figure 7a is the log likelihood functionL(θ) of the first 8 natural frequencies when the parameter settings of the test data samples are taken as σ E = 5 GPa, L E = 500 mm. It can be seen that, whenL(θ) attains its maximum value, the parameter θ * of the corresponding point is also (5 GPa, 500 mm), which shows that the random field model introduced can identify the model parameter very well, and the random field model is reliable.
In Figure 7b, based on Equation (32), a group of parameter setting values of test samples are randomly taken as σ E = 6.8 GPa, L E = 210 mm, and the multi-dimensional kernel density probability distribution function is estimated for the first 8 natural frequencies, and then the maximum likelihood estimation is carried out. The estimated parameter results, that is θ * = (7 GPa, 200 mm), are slightly different from the original parameter setting (6.8 GPa, 210 mm). This is because the amount of the model samples may not be infinite and the estimation accuracy of the input parameter of test samples is obviously limited by the amount of the model samples, but the peak value of likelihood function can still be obtained around the input parameter θ = (6.8, 210).
When considering random fields E and ρ simultaneously, the input parameters of the test samples are taken as two groups, respectively: θ E = (σ E , L E ) = (5 GPa, 500 mm), θ ρ = (σ ρ , L ρ ) = (250 kg/m 3 , 500 mm), as well as θ E = (σ E , L E ) = (4 GPa, 500 mm), θ ρ = (σ ρ , L ρ ) = (450 kg/m 3 , 500 mm). The correlation length of two random fields is fixed; that is, L E = L ρ = 500 mm; the mean variances σ E and σ ρ are then estimated when the likelihood functionL(θ) is taken as its maximal value. It can be seen from Figure 8 that the parameter setting at the maximum value point ofL(θ) is the same as the input parameters by using the maximum likelihood estimation for the first 8-order natural frequency.

The Expansion and Distribution of Two-Dimensional Random Field
First, considering E is a two-dimensional random field, E is then expanded with K-L expansion as ( ) ( )  Figure 10.

The Expansion and Distribution of Two-Dimensional Random Field
First, considering E is a two-dimensional random field, E is then expande K-L expansion as ( ) ( )  Figure 10.

The Expansion and Distribution of Two-Dimensional Random Field
First, considering E is a two-dimensional random field, E is then expanded with K-L expansion as Taking the parameter setting of random field E, θ E = σ E , L Ex , L Ey = (5GPa, 500mm, 500mm), as the input parameter of test samples and model samples. When the plate is meshed into different amounts of elements, the distribution of random field E in the plate is displayed in Figure 10.
Moreover, only considering mass density ρ as random field, ρ is expressed with the two-dimensional (2D) K-L expansion as Equation (37), whereby its mean value is ρ(x, y) = 7800kg/m 3 . Taking its parameter setting as θ ρ = σ ρ , L ρx , L ρy = (7800kg/m 3 , 300mm, 300mm) and meshing the plate into different number of elements, the distribution of ρ in the steel plate is displayed in Figure 11. Moreover, only considering mass density ρ as random field, ρ is expresse the two-dimensional (2D) K-L expansion as Equation (37), whereby its mean v

Kernel Density Estimation and Maximum Likelihood Estimation
When only considering random field E , a group of parameter settings are ta   Moreover, only considering mass density ρ as random field, ρ is expressed with the two-dimensional (2D) K-L expansion as Equation (37), whereby its mean value is

Kernel Density Estimation and Maximum Likelihood Estimation
When only considering random field E , a group of parameter settings are taken as It can be seen from Figures 10 and 11 that the values of E and ρ fluctuate and vary around the mean values E(x, y) and ρ(x, y) , and the values of E and ρ randomly distributed along x and y directions. Moreover, compared with the random distributions of E and ρ plated with 225 meshing elements, the random distributions with 400 elements obviously reflect the real cases better and more accurately.
From Table 1, according to the point at whichL(θ) attains its maximal value, the parameter setting of test samples, i.e., θ * E = σ * E , L * Ex , L * Ey , can be well estimated by the presented 2D random field model of the plate, which is basically consistent with the input parameters θ E1 , θ E2 , θ E3 , and θ E4 of the test samples. Hence, the constructed 2D random field model is reliable. Similar to the 1D random field model of I-beam, the amount of model sample groups may not be infinite, and the estimation accuracy of test sample parameter is limited by the amount of the model samples. When the maximum likelihood estimation method is used to estimate the test samples with input parameter θ E4 = (9GPa, 430mm, 650mm), the input parameter of test samples can still be estimated comparatively accurately; i.e., θ * E4 = (9GPa, 400mm, 700mm). Furthermore, when considering random fields E and ρ simultaneously, their parameter settings are, respectively, taken as θ E = σ E , L Ex , L Ey = (σ E , 500mm, 500mm) and θ ρ = σ ρ , L ρx , L ρy = (250kg/m 3 , 500mm, 500mm). When σ E is taken as 1GPa, 2GPa, . . . , 10GPa, respectively, parameter settings of the test samples are taken as θ E0 = (5GPa, 500mm, 500mm) and θ ρ0 = 250kg/m 3 , 500mm, 500mm ; the 2D kernel density distribution function of the second-order natural frequency f 2 of test samples is estimated as shown in the left part of Figure 12. Similarly, parameter settings of E and ρ are θ E = σ E , L Ex , L Ey = (5 GPa, 500 mm, 500 mm) and θ ρ = σ ρ , L ρx , L ρy = σ ρ , 500 mm, 500 mm , and, respectively, taking σ ρ as 50 kg/m 3 , 100 kg/m 3 , . . . , 500 kg/m 3 for computing, and then taking parameter settings of test samples θ E0 = (5 GPa, 500 mm, 500 mm) and θ ρ0 = 250 kg/m 3 , 500mm, 500 mm , the 2D kernel density distribution function of f 2 is estimated as shown in the right part of Figure 12.
When taking input parameter settings of test samples θ E0 = (5 GPa, 500 mm, 500 mm) and θ ρ0 = (250 kg/m 3 , 500 mm, 500 mm), the curve of kernel density estimation of f 2 and the histograms of probability density of f 2 are illustrated in Figure 13.   Figure 12 shows again that random field E has a greater influence on the kernel density distribution function of structural natural frequency than ρ does. In Figure 13, the curve of kernel density estimations and the histograms of 2 f agree well, and it is obvious that the curve and histogram in the right part are more reasonable than in the left part with the increasing amount of meshing elements.  Figure 12 shows again that random field E has a greater influence on the kernel density distribution function of structural natural frequency than ρ does. In Figure 13, the curve of kernel density estimations and the histograms of f 2 agree well, and it is obvious that the curve and histogram in the right part are more reasonable than in the left part with the increasing amount of meshing elements.
In addition, fixing mean variance of parameter θ E , σ E = 5GPa, and, taking L Ex and L Ey as variables, the variation in log likelihood functionL(θ) of natural frequency with L Ex and L Ey is shown in Figure 14a. Similarly,L(θ) is obtained in Figure 14b when parameter settings of test samples σ E = 2.1GPa, L E = L Ex = L Ey = 300mm. In Figure 14, the plate is meshed into 400 elements.  Figure 14a. Similarly,  ( ) L θ is obtained in Figure 14b when parameter settings of test samples In Figure   14, the plate is meshed into 400 elements.
(a) 500mm 5GPa, It can be seen from Figure 14 that the parameter settings  It can be seen from Figure 14 that the parameter settings θ * E obtained corresponding to the maximal valueL(θ) are, respectively, θ * E = (σ * e , L * Ex , L * Ey ) = (5GPa, 500mm, 500mm) and (2GPa, 300mm, 300mm); based on the multi-dimensional kernel density estimation, parameter settings of test samples can be accurately estimated when the log likelihood function attains its maximal value and the estimated parameter settings are very close to the input parameters of the test samples, which verifies the validity of the constructed model.
In the same way, when only ρ is random field and the input parameter settings θ ρ of the test samples are σ ρ = 500 kg/m 3 , L ρx = 400 mm, L ρy = 700 mm and σ ρ = 150 kg/m 3 , L = L ρx = L ρy = 600 mm, respectively,L(θ) obtained based on the first six natural frequencies of the plate is displayed in Figure 15. 500mm , the distributions of the lower bound, mean value, and upper bound of the first two random natural modes are illustrated in Figure 17.  Figure 16 shows the distributions of the lower bound, mean value, and upper bound of the first two-order random natural modes in the steel plate when only considering random field E, and the parameter setting after K-L expansion is taken as θ E = σ E , L Ex , L Ey = (5 GPa, 500 mm, 500 mm). When considering random fields E and ρ simultaneously and taking parameter settings θ E = (σ E , L Ex , L Ey ) = (5GPa, 500mm, 500mm) and θ σ = σ σ , L ρx , L ρy = (250kg/m 3 , 500mm, 500mm), the distributions of the lower bound, mean value, and upper bound of the first two random natural modes are illustrated in Figure 17.  Table 2 for different random models so as to compare the influences of different random cases on the structural random natural modes.

Investigation of the Random Characteristics
(a) The first order natural mode (b) The second order natural mode of random natural modes when considering random fields E and ρ simultaneously. Table 2. Comparison of the norm of natural modes for different random models.  Table 2, it can be seen that mean values of random natural modes only considering random field E are very close to those simultaneously considering random fields E and ρ , but mean variances for these two random models are very different and mean variances only considering the randomness of E are obviously smaller than those considering the randomness of E and ρ simultaneously.

Norm of Random Natural Modes
The mean values  Table 3 for different random models so as to compare the influences of different random cases on the structural random natural frequencies.
Similarly, investigation of the influences of different random models on random natural frequencies is implemented and the corresponding results are illustrated in Figure  18. From Figure 18 and Table 3, once again, the mean variances and value ranges of random natural frequencies simultaneously considering the randomness of E and ρ are obviously greater than those only considering the randomness of E , and the former reflects the more realistic case in structural engineering than the latter does. The mean values Φ i and mean variances σ Φ i of norms of the first 4-order natural modes are computed and listed in Table 2 for different random models so as to compare the influences of different random cases on the structural random natural modes. Table 2. Comparison of the norm of natural modes for different random models.  Table 2, it can be seen that mean values of random natural modes only considering random field E are very close to those simultaneously considering random fields E and ρ, but mean variances for these two random models are very different and mean variances only considering the randomness of E are obviously smaller than those considering the randomness of E and ρ simultaneously.

Norm of Random Natural Modes
The mean values f i and mean variances σ f i of the first four natural frequencies are computed and listed in Table 3 for different random models so as to compare the influences of different random cases on the structural random natural frequencies. Similarly, investigation of the influences of different random models on random natural frequencies is implemented and the corresponding results are illustrated in Figure 18.
From Figure 18 and Table 3, once again, the mean variances and value ranges of random natural frequencies simultaneously considering the randomness of E and ρ are obviously greater than those only considering the randomness of E, and the former reflects the more realistic case in structural engineering than the latter does.

Conclusions
In this work, an investigation on the stochastic free vibration problem of eng structures considering material uncertainties is presented. As a novel extension o ventional uncertain eigenvalue problem, spatially dependent stochastic parame random field theory are combined into a numerical analysis framework of stocha element method, and the verification method to validate the proposed paramete model and stochastic free vibration model is presented and updated by using t mum likelihood method: (1) The parameter setting model based on random field theory can represent tially dependent uncertainty of structural parameters well, and the paramet model presented can describe the randomly varying characteristics of actu tural parameters. (2) The example shows that the parameter settings of the model can be quantifi output response of the structural system; i.e., structural dynamic characteris as the structural natural frequency, and the mean variance and autocorrela tance of the parameter of the structure can also be obtained, which is very im to application of random field in engineering. Figure 18. The value range of random natural frequencies for different random models computed by f n ± 3 × σ f n .

Conclusions
In this work, an investigation on the stochastic free vibration problem of engineering structures considering material uncertainties is presented. As a novel extension of the conventional uncertain eigenvalue problem, spatially dependent stochastic parameters and random field theory are combined into a numerical analysis framework of stochastic finite element method, and the verification method to validate the proposed parameter setting model and stochastic free vibration model is presented and updated by using the maximum likelihood method: (1) The parameter setting model based on random field theory can represent the spatially dependent uncertainty of structural parameters well, and the parameter setting model presented can describe the randomly varying characteristics of actual structural parameters.
(2) The example shows that the parameter settings of the model can be quantified by the output response of the structural system; i.e., structural dynamic characteristics, such as the structural natural frequency, and the mean variance and autocorrelation distance of the parameter of the structure can also be obtained, which is very important to application of random field in engineering. (3) The proposed method can be extended to apply to other structural parameters and can also be used to establish and quantify the parameter setting model of random fields for other material parameters or structural parameters. The applicability and effectiveness of the proposed computational framework are evidently demonstrated through the numerical investigations on various practically motivated engineering structures. (4) Obviously, the simulation results are closer to reality when more parameters are considered with the random field. However, as the number of parameters considered increases, the computational effort increases exponentially. How to strike a valuable trade-off between them is an interesting area of future work.