An Interpolation Method for Karhunen–Loève Expansion of Random Field Discretization

Zi Hanand Zhentian Huang

1Department of Engineering Mechanics,Hohai University,Nanjing,211100,China

2College of Civil and Transportation Engineering,Shenzhen University,Guangdong Provincial Key Laboratory of Deep Earth Sciences and Geothermal Energy Exploitation and Utilization,Institute of Deep Earth Sciences and Green Energy,Shenzhen,China

ABSTRACT

In the context of global mean square error concerning the number of random variables in the representation,the Karhunen–Loève(KL)expansion is the optimal series expansion method for random field discretization.The computational efficiency and accuracy of the KL expansion are contingent upon the accurate resolution of the Fredholm integral eigenvalue problem (IEVP).The paper proposes an interpolation method based on different interpolation basis functions such as moving least squares(MLS),least squares(LS),and finite element method(FEM)to solve the IEVP.Compared with the Galerkin method based on finite element or Legendre polynomials,the main advantage of the interpolation method is that,in the calculation of eigenvalues and eigenfunctions in one-dimensional random fields,the integral matrix containing covariance function only requires a single integral,which is less than a two-folded integral by the Galerkin method.The effectiveness and computational efficiency of the proposed interpolation method are verified through various one-dimensional examples.Furthermore,based on the KL expansion and polynomial chaos expansion,the stochastic analysis of two-dimensional regular and irregular domains is conducted,and the basis function of the extended finite element method(XFEM)is introduced as the interpolation basis function in two-dimensional irregular domains to solve the IEVP.

KEYWORDS

Random field discretization;KL expansion;IEVP;MLS;FEM;stochastic analysis

1 Introduction

In civil engineering,it is commonplace for the material properties and geometric dimensions of structures to exhibit inherent uncertainties.In the context of engineering systems,it is imperative to consider the presence of uncertainties during the analysis and design phase.These uncertainties are commonly attributed to the random spatial variation of certain properties.Traditional approaches that rely on simulating the characteristics of random structures with single random variables are deemed insufficient,as they fail to account for point-to-point variability [1].An alternative mathematical model that has gained significant traction in the field is the continuous random field model,which accurately captures the spatial variation of a structure’s properties.Specifically,a random field is a stochastic process defined on a continuous region,where each point is associated with a random variable.However,the computation of the entire random field is unfeasible due to the infinite number of random variables involved.To facilitate the calculations,a limited number of random variables must be applied to represent the random field,which is referred to as random field discretization[2].

The methods used for random field discretization can be broadly categorized into two types:spatial discretization and series expansion.Spatial discretization methods comprise the central point method[3],the local mean method[4],the shape function method[5],and the optimal linear estimation method [6].Among these,the local average method is commonly used due to its stable calculation results and insensitivity to the structures associated with the random field.The efficiency and accuracy of the spatial discretization methods are closely related to the random field mesh.

Series expansion methods used for random field discretization include Karhunen–Loève (KL)expansion [7,8],orthogonal polynomial expansion [9],etc.The KL expansion,introduced by Ghanem et al.[10],has emerged as the most extensively applied series expansion method for random field discretization,particularly in the context of uncertainties related to input parameters [11].Moreover,when considering the global mean square error with respect to the number of random variables,the KL expansion is optimal compared to other series expansion methods [10].KL expansion uses a linear combination of orthogonal bases to represent random fields,and the selected orthogonal functions are the eigenfunctions of the Fredholm integral equation of the second kind.The integral kernel function in the Fredholm integral equation is the covariance function of the random field.For simple geometries and special covariance functions,analytical solutions can be obtained.However,for more general scenarios,numerical methods are required to solve a Fredholm integral eigenvalue problem(IEVP),with the Galerkin method being the most commonly used[2].

In one-dimensional domains,the Galerkin method is frequently employed in a spectral sense with basis functions spanning over the entire domain.The convergence behavior of this approach is examined in[12]for both stationary and non-stationary problems using polynomials with order up to ten as basis functions.Other polynomial bases,such as Legendre polynomials[9],Gegenbauer polynomials[13],and Chebyshev polynomials[14],have also been proposed.However,these methods can only improve the calculation accuracy by increasing the polynomial order,which requires significant computational effort.The wavelet-Galerkin method[15]also represents a Galerkin-based technique for solving Fredholm integral eigenvalue problems in one-dimensional domains.For two-dimensional and three-dimensional domains with random fields,Ghanem et al.[10] advocated using the finite element method (FEM) to obtain approximate solutions for the KL expansion,while Papaioannou[13] examined the convergence of the FEM in two-dimensional domains.Based FEM-Galerkin method for the discretization of the IEVP,Allaix et al.[1] proposed a genetic algorithm to achieve an optimal discretization of 2D homogeneous random fields.To compute matrix eigenvalue problems in the FEM,a generalized fast multi-pole Krylov eigen-solver is employed [16].In two-dimensional and particularly three-dimensional random fields,the computational cost of constructing the matrix eigenvalue problem and computing its solution can be prohibitively expensive.Basmaji et al.[17]proposed a Galerkin scheme based on discontinuous Legendre polynomials to solve IEVP,which reduces the computational complexity by constructing the Legendre basis on each local element domain and exploiting the orthogonality of Legendre polynomials.In addition,generating meshes for complex geometries may also present difficulties and time constraints.To circumvent this issue,Rahman et al.[18]employed meshless basis functions in the Galerkin method for complex domains.Papaioannou [13] introduced a spectral Galerkin method that addresses the challenges posed by geometrically complex domains by embedding the actual domain into a simple geometric shape.Rahman [19] proposed the Galerkin isogeometric method for KL approximation which can denote computational domains by non-uniform rational B-splines(NURBS).

The development of a method that can efficiently solve the IEVP with reduced computational requirements and without the need for a complex mesh is crucial for effectively addressing complex domain problems.In this study,an interpolation approach for solving the IEVP in KL expansion is proposed.And the impact of three different interpolation basis functions,namely moving least squares [20–22] (MLS),least squares [23] (LS),and finite element method (FEM),on the accuracy of computation is explored.MLS and LS are widely used in meshless methods[22,24–26].The results demonstrate that for one-dimensional domains,the proposed interpolation method,which requires a single integral for the integral matrix containing the covariance function,is more computationally efficient than the Galerkin method,which requires a two-folded integral.

In addition,we employ a combination of KL expansion and polynomial chaos expansion to perform stochastic analysis in two-dimensional domains.To tackle the issue of grid division in irregular domains,we introduce the extended finite element method(XFEM)as an alternative to the traditional finite element method.The XFEM[27,28]is a well-known numerical method that addresses discontinuity problems based on the partition of unity theory [29].This approach adds enrichment functions to the displacement mode of conventional finite element methods to accurately represent discontinuities.The key advantage of the XFEM is that the computational mesh is independent of the internal geometry or physical interface of the structure,making it an effective method for solving problems that involve discontinuities,especially in the presence of holes[30].

The structure of this paper is as follows: Section 2 introduces the KL expansion and its fundamental properties.Additionally,the numerical methods used to solve IEVP,namely the Galerkin method and the interpolation method,are discussed.Interpolation basis functions such as MLS,LS,and FEM are provided.Section 3 briefly describes the spectral stochastic finite element method and extended finite element method,which are then used for stochastic analysis of complex domains.Section 4 is dedicated to numerical examples,wherein the computational efficiency of the numerical methods is compared in one-dimensional domains,and the stochastic analysis is conducted in both 2D regular and irregular domains.The benefits of the proposed method are highlighted.At last,Section 5 concludes the paper.

2 Karhunen–Loève Expansion

Defining(Θ,F,P)is a complete probability space,Θis the sample space.A continuous random fieldH(x,θ)is a measurable functionH:Ω×Θ→Rindexed by the space coordinatesx∈Ω∈Rn,whereΩis a continuously bounded domain.For a givenx0∈Ω,H(x0,θ) is a random variable.For a given outcomeθ0∈Θ,H(x,θ0) is a realization of the field.In practice,it is difficult to apply a continuous random field directly,so the random field needs to be treated with discretization.Karhunen–Loève expansion is a series expansion method for representing the random field.which is founded upon the spectral decomposition of the field’s covariance function.The Karhunen–Loève expansion of a second-order random field can be expressed as[10]

whereμ(x) is the mean function,ξi(θ) are standard uncorrelated random variables,i.e.,E[ξi(θ)]=0,E[ξi(θ)ξj(θ)]=δij.λiandφiare the eigenvalues and eigenfunctions of the covariance function respectively which can be obtained from solving the homogeneous Fredholm integral equation of the second kind:

where Cov(x,x′)is the covariance function which can be represented as Cov(x,x′)=σ(x)·σ(x′)·ρ(x,x′),σis the standard deviation function of the random field andρis the correlation function.Cov(x,x′) is also called as kernel function,and it is bounded,symmetric and positive semi-definite[4].

According to Mercer’s Theorem [31],the eigen-decomposition of the covariance function is as follows:

The eigenvaluesλiare nonnegative,and the eigenfunctionsφi(x) form a completely orthogonal set satisfying the following equation:

whereδijis the Kronecker-delta function.

Due to the orthogonality of eigenfunctions,it is easy to get a closed form of each random variable as follows:

In the case where the random fieldH(x,θ)is Gaussian,theξiare independent standard normal random variables [10].The joint distribution ofξiis almost impossible to acquire in other cases.However,scholars have also studied the non-Gaussian random process simulation based on KL expansion.Phoon et al.[32] developed an iterative framework to simulate non-stationary and non-Gaussian processes.Sakamoto et al.[33]used the polynomial chaos expansion of a Gaussian stochastic process to represent the non-Gaussian process,while the Gaussian stochastic process is modeled by the KL expansion.Dai et al.[34]used the one-dimensional polynomial chaos expansion to decompose the random coefficients in the KL expansion and obtained the explicit representation of non-Gaussian and non-stationary random processes described by the covariance function and the marginal cumulative distribution function which is suitable for random finite element analysis of structures.Tong et al.[35]combined Karhunen–Loève expansion with the Linear-moments-based(L-moments-based)Hermite polynomial model for simulating strongly non-Gaussian and non-stationary processes.

For practical implementation,the random fieldH(x,θ) can be acquired through sorting the eigenvaluesλiin a descending order and truncating the series expansion at theM-th term

The second-order statistics of the truncated series Eq.(6)can be presented as

AsM→∞,Eq.(7) converges to Eq.(3).The variance function of the approximated random field is

Integration of the expectation of the squared approximation error over the domainΩyields the global mean square error[10]

An alternative error estimation is the error variance[6]which is defined as follows:

The corresponding mean error variance is given as

where|Ω|=

If the variance of the random field is constant,i.e.,∀x∈Ω,σ(x)=σ,the mean error variancecan be reduced to

The error variance in this section is based on the analytical solutions of eigenvalues and eigenfunctions.However,it is not always possible to obtain the analytical eigen-solutions in practical engineering applications.Therefore,in the following sections,we will introduce other error measures in the numerical examples to evaluate the accuracy of the proposed interpolation method.

2.1 Numerical Solution of the Fredholm Integral Equation

The key to KL expansion is to find the eigenvalues and eigenfunctions of the covariance function by solving a Fredholm integral eigenvalue problem of Eq.(2).The analytic solution of Eq.(2) is only applicable to a finite set of covariance functions and geometric domains.For general cases,the numerical method is the only recourse.Therefore,the random field approximation of truncated KL expansion in Eq.(6)can be approximated as

2.1.1TheGalerkinMethodforFredholmIntegralEquation

Let {hk(x)} be a complete basis of the Hilbert spaceL2(Ω).In the Galerkin method,the eigenfunction(x)can be reprsented over the basis as

whereare the coefficients of thei-th eigenfunction.

Substituting Eq.(15) into Eq.(2) and adopting a Galerkin procedure for the corresponding residual,the equation is obtained as follows:

Eq.(16)could rewritten as

where

The eigenvaluesλiexists in the diagonal of matrixΛ.By solving the generalized matrix eigenvalue problem of Eq.(17),the approximate eigenvalueand the coefficients of basis functiondikare obtained.Substitutingdkito Eq.(15)can get approximate eigenfunctions(x).

Various types of basis functions can be utilized with the Galerkin method,including orthogonal polynomials such as Legendre polynomials,which yield a diagonal matrixG.However,increasing the order of these polynomials is the only way to enhance the accuracy of the approximation,and this can lead to numerical issues if the correlation length of the random field is too small.Another frequently used basis function is the shape function of the finite element method,also known as the FEM-based Galerkin method[2].However,this method requires domainΩmeshing,which can be a challenging task for complex geometries.Moreover,the computation of the matrixCcan be time-consuming,particularly in two and three dimensions.

2.1.2TheInterpolationMethodforFredholmIntegralEquation

In the random field discretization based on KL expansion,it is necessary to integrate all the elements in matrixCby using the Galerkin method to solve IEVP.In the case of a one-dimensional domain,the elements inCinvolve a double integral,whereas for a two-dimensional domain,a fourfolded integral is required.As the dimensionality of the problem increases,the computational cost escalates significantly.To address this issue,we propose an interpolation method for solving the IEVP of Eq.(2)which is a homogeneous Fredholm integral equation of the second kind.To better illustrate the computing process of the proposed interpolation method,we start with the general form of the Fredholm integral equation of the second kind,

whereP(x),Q(x),f(x) are continuous functions on the interval [a,b],andP(x) /=0.The integral kernel functionL(x,x′)is a continuous function of two variables on the interval[a,b]×[a,b].

The key problem of the Fredholm integral equation is the discretization of unknown functions.Firstly,the integral interval[a,b]is discretized asa=x1＜x2＜···＜xn=b,n is the number of nodes.m1,m2,...,mnare the corresponding function values of nodes.Then the approximating functionm(x)can be expressed as

whereNI(x)is the basis function of interpolation,andIis the node of interpolation.

Substituting Eq.(23)into Eq.(22),and let the equation be true at the nodesxJ(J=1,2,...,n),one can obtain

Switching the order of integration and summating the integral term of right-hand side of Eq.(24)yield

where

Eq.(24)is rewritten as

IfI=J,NI(xJ)=1,I/=J,NI(xJ)=0,Eq.(27)can be presented in a matrix form as

where

Comparing Eqs.(2) and (22),we can see thatf(x)=0,Qis a identity matrix andP=λI.Therefore,Eq.(28)can be converted into the following form:

By solving the eigenvalue problem of Eq.(33),the approximate eigenvaluesand the corresponding node valuesmcan be calculated.The eigenfunctions can be obtained by substitutingminto Eq.(23).

Compared with the Galerkin method,which requires a two-folded integral to calculate the elements of matrixCin one dimension,the interpolation method only requires a single integral to calculate the elements of matrixL.Obviously,interpolation method needs less computation.

The selection of the basis function is a crucial step in interpolation methods.This paper discusses three choices of basis functions,namely moving least squares[20–22](MLS),least squares[23](LS),and finite element method(FEM).The MLS is primarily introduced,and the LS can be derived from it.Although the finite element method is not discussed in detail in this paper,it is another viable option for basis functions.

If the domain Ω discretes intonnodes and the values of eigenfunction at the nodesxIare known.The approximate eigenfunction ˆφ(x)using moving least squares interpolation can be written as follows:

wheretis the number of terms in the basis,p(x)is a polynomial basis,a(x)is the vector of coefficients,which are functions of the spatial coordinatesx.For 1D problems,commonly used bases arep=[1,x]Tfort=2 andp=[1,x,x2]Tfort=3,etc.In the same way,for 2D problems,generally used bases arep=[1,x,y]Tfort=3 andp=[1,x,y,x2,xy,y2]Tfort=6,etc.

The weighted square sum of error of functionm(x)at nodes is

wheremI=m(xI)are the values of function at nodexI,wI(x)is a weight function that is greater than zero only in a finite domiam ΩIaround nodex,and it is zero outside of ΩI,as shown in Fig.1.wI(x)is compactly supported and ΩIis the support domain.A typical choice forwI(x) is the normalized Gaussian function,namely

whereβis a constant,r=‖x-xI‖/dmI,anddmI=ds×cIis the radius of the support domain of pointxI,dsis the multiplier greater than 1,andcIis the distance between the nodexIand its nearest node.

Minimizing the weighted square sum of errorJ,namely

Figure 1 :The support domains of nodes

Then it leads to

where

The vectora(x)can be obtained from Eq.(38),and substituting it into Eq.(34).Then the shape function of moving least squares can be obtained as

In the process,if the weight functionωI(x)is equal to

The moving least squares is reduced to least squares(LS),and matrixAof Eq.(39)and matrixBof Eq.(40)become

3 The Stochastic Analysis

In this section,we first describe the theoretical formulation of the spectral stochastic finite element method which is a numerical approach to model the random parameter system in terms of finite element framework.While,for the irregular domains,we applied the extended finite element method instead of the finite element method for calculation.It is well-known that the extended finite element method(XFEM)is an effective numerical method for solving discontinuity problems[27,28].Based on the partition of unity method[29],XFEM adds some enrichment functions into the displacement mode of the finite element method to reflect discontinuity,and it has the advantage that the mesh is independent of the geometry or physical interface inside the structure.Therefore,for random analysis of irregular domains,XFEM has more advantages.Moreover,the basis function of XFEM can also be used as the basis function of the interpolation method to solve a Fredholm integral eigenvalue problem in KL expansion.

The spectral stochastic finite element method introduced by Ghanem et al.[10]is an extension of the deterministic finite element method for solving boundary value problems of stochastic material properties.It supposes that the material Young’s modulus is a Gaussian random field.The elasticity matrix in pointxis written as

whereD0represents a constant matrix.

Applying the KL expansion in Eq.(1),the stochastic matrix of a finite element has the following form:

whereke0is the mean element stiffness matrix,keiare deterministic matrices given by

whereBis the strain-displacement matrix.

Assuming that the loading is determined and denotingξ0(θ)=1,the finite element equilibrium equation is

where thee nodal displacement vectoruis presented by polynomial chaos expansion(PCE)as follows:

Substituting the expansion Eq.(50)into Eq.(49)yields:

For calculating purposes,truncating the KL expansion afterMterms and polynomial chaos expansion afterPterms results in

whereP=(p+M)!/(p!M!),pis the order of polynomial chaos.

By making the residual orthogonal to the approximate space spanned by PCE,one can get

wherecijk=,andFk=E(ΨkF).E(·)denotes the mathematical expectation.

Eq.(53)can be rewritten as

where

Eq.(54)is a system of equations of sizeNP×NP,Nis the number of physical degrees of freedom.EachKjkis a matrix of sizeN×N.Eachujis a N-dimensional vector.As the coefficient vectorsujis obtained,the mean and covariance matrix ofu(θ)can be obtained as

Substituting the displacement form of the finite element with that of the XFEM can solve the elastic modulus stochastic problem in the irregular domain.In XFEM,the displacement approximation is expressed as[27,36]

whereNiis the shape function as used in the formulation of conventional FEM;ui,ajare the nodal displacements and nodal enrichment variables,respectively;Φ(x)is the enrichment function associated with the discontinuity;ISis the set of all nodes in the discrete mesh;IEis the set of nodes of the elements which contain the interface.

For models with holes,the XFEM displacement approximation has a simple form as below[30]:

whereHis the Heaviside function and the value is as follows:

whereψ(x)is the level set.In general,the signed distance function is adopted to represent the level set function,and is defined as

wherexiis thei-th node coordinates,and ‖ · ‖ is theL2norm,xΓis the interface of the hole.If the hole is circular,the level set can be represented asψ(x)=‖x-c‖-R,xis the coordinate of the node in the extended finite element mesh,c=(xc,yc)is the center of the hole,‖ ·‖represents the distance between the node and the center of the hole,andRis the radius of the hole.

4 Numerical Studies

4.1 Error Measures

In the case of random fields without analytical eigen-solutions,the error estimation methods outlined in Section 2 cannot produce computable expressions.As a result,two error measures are introduced[14]to enable numerical evaluation of the interpolation method and Galerkin method.

If the Fredholm integral eigenvalue problem of KL expansion has analytical solutions,(x,θ)in Eqs.(62)and(63)can be replaced by(x,θ).

4.2 Discretization of One-Dimensional Random Fields

To demonstrate the validity and advantages of the interpolation approach for solving Fredholm integral eigenvalue problems in KL expansion,the method is applied to a variety of covariance functions and compared against the conventional Galerkin method.

Example 1 is a homogeneous random field,with the exponential kernel function defined as

wherex∈[-a,a],x′∈[-a,a],bis the correlation length,σis the standard deviation of the random field,andσ=1,a=1,b=1.

Example 2 is a non-homogeneous random field,with the Wiener-Lévy kernel function defined as

wherex∈[0,a],x′∈[0,a],a=1 andσ=1.

Example 3 is a random field with squared exponential covariance function defined as

where the correlation lengthb=1,x∈[-a,a],a=1 and the standard deviation of the random field isσ=1.

The eigenvalues and eigenfunctions corresponding to the kernel functions in examples 1 and 2 have analytical solutions[37].The accuracy of different numerical methods for solving Fredholm integral eigenvalue problems in KL expansion was compared with that of the analytical solutions.The error of eigenfunctions is defined as follows:

In solving Fredholm integral eigenvalue problems using the interpolation method,we select interpolation basis functions including moving least squares (MLS),least squares (LS),and finite element method(FEM).In Fig.2 displaying computed results,these functions are denoted as MLS,LS,and FEM-interpolation,respectively.For comparison purposes,the finite element basis function and the Legendre polynomial are applied as the basis functions of the Galerkin method and denoted as FEM-Galerkin and Legendre,respectively,in the figures.The number of nodes in the computation of examples 1 and 2 is set toNn=100.In KL expansion,the number of truncations isM=20.In MLS,the scale parameter is set tods=1.5,and the basis function ist=2.In LS,the scale parameter is set tods=0.9.The order of the basis function in FEM is 1,and the order of the Legendre polynomial isp=21.

Figure 2 :The eigenvalues(a)and the errors of eigenfunctions(b)of the exponential kernel function

Figs.2a and 3a display the eigenvalues of examples 1 and 2,respectively,while Figs.2b and 3b show the errors of eigenfunctions of examples 1 and 2,respectively.The computational formula is presented in Eq.(67).Fig.3 indicates that the eigenvalues calculated using the interpolation method based on MLS,LS,and FEM are consistent with the analytical solutions.Additionally,the accuracy is close to that of the Galerkin method based on FEM,while the accuracy of the Galerkin method based on Legendre polynomials is slightly poor.Concerning the eigenfunctions,the accuracy of the FEMbased interpolation method is superior to or close to other methods,and it hardly fluctuates with an increase in the eigenfunction index.The accuracy of eigenfunctions calculated by the interpolation method based on MLS and LS is similar.In example 1,the error of eigenfunctions calculated using the Galerkin method based on the Legendre polynomial decreases rapidly as the eigenfunction index increases.

To investigate the impact of various correlation lengths on eigenvalues,different correlation lengths(b=0.4,1,1.6,etc.)were selected for examples 1 and 3,and the results are presented in Fig.4.The comparison of results presents that the correlation length has a more significant influence on the first few eigenvalues,and the longer the correlation length,the steeper the slope of the first few terms of the eigenvalue curve.The computed results of the MLS,LS,and FEM-interpolation are similar,indicating that the correlation length does not influence the accuracy of eigenvalues calculated using the interpolation method.

Figure 3 :The eigenvalues(a)and the errors of eigenfunctions(b)of the Winer-Levy kernel function

Figure 4 : The eigenvalues of the exponential kernel function (a) and squared exponential kernel function(b)under different correlation lengths b

The preceding discussion confirms the accuracy of the methods.The subsequent discussion focuses on examining the impact of various parameters in the interpolation method on the computed results.For the FEM-based interpolation method,the relationship between the computed results and the number of nodes is investigated.Fig.5 illustrates the mean of relative variance error∈Var(a)calculated from Eq.(62)and the mean of relative covariance error∈Cov(b)calculated from Eq.(63)for examples 1 and 2,respectively.The computed errors are compared with∈Varand∈Covcalculated from the analytical eigenvalues and eigenfunctions.The results demonstrate that∈Varand∈Covcalculated using the FEM-based interpolation method gradually converge as the number of nodes increases,when the number of nodesNnis close to 120,the errors tend to converge,and∈Covis less than∈Var.

In the case of MLS-and LS-based interpolation methods,this study mainly discusses the influence of the number of nodesNnand the multiplierdsin the node support domain on the computed results.Figs.6 and 7 show the mean of relative variance error∈Var(a)and the mean of relative covariance error∈Cov(b)calculated using the MLS-based interpolation method in examples 1 and 2,respectively.The results indicate that as the number of nodesNnincreases,both∈Varand∈Covgradually converge,and∈Covis less than∈Var.When the number of nodes is small,the multiplierdshas a significant influence on the error.However,as the number of nodes increases,the influence of the multiplierdsdecreases rapidly,when the number of nodesNnis close to 120,the influence of the multiplierdstends to be small.

Fig.8 displays the mean of the relative variance error∈Var(a) and the mean of the relative covariance error∈Cov(b)calculated using the MLS-based interpolation method in example 3.In this example,the influence of the multiplierdsof the node support domain on the errors significantly increases.The smaller the multiplierds,the lower the mean of relative variance error∈Var(a)and the mean of relative covariance error∈Cov(b),when the multiplierdsof the node support domain is closer to 1,the computational accuracy is higher.Under differentds,along with the increase of the number of nodes,∈Varand∈Covgradually decrease,the values of∈Varand∈Covin this example are two orders of magnitude smaller than those in examples 1 and 2,whenNn=120,ds=1.01,the values of∈Varand∈Covcan reach 10-4.

Figure 7 :The convergence curves of the MLS-based interpolation method with respect to the number of nodes Nn in example 2,the mean of relative variance error ∈Var (a),the mean of relative covariance error ∈Cov(b)

Figure 8 :The convergence curves of the MLS-based interpolation method with respect to the number of nodes Nn in example 3,the mean of relative variance error ∈Var(a),the mean of relative covariance error ∈Cov(b)

Figs.9–11 display the mean of relative variance error∈Var(a)and the mean of relative covariance error∈Cov(b)calculated using the LS-based interpolation method in examples 1,2 and 3,respectively.It is evident that when the multiplierdsof the node support domain is small,increasing the number of nodes does not enable the error to gradually converge.However,when the multiplierdsapproaches 1,∈Varand∈Covgradually converge with the increase of the number of nodesNnexcept for the individual points of example 3.WhenNn=120,ds=0.96,∈Varand∈Covtend to converge or obtained relatively good accuracy.Based on the comparison of the MLS-and LS-based interpolation methods,it is known that,in some examples,the MLS-based interpolation method is less sensitive to the multiplierdsof the node support domain,while the LS-based interpolation method is relatively sensitive and the multiplierdsshould as close to 1 as possible.The reason is that,in the LS-based interpolation method,the moredsapproaches 1,the moreNI(xJ)approaches 1 whenI=J.It is consistent with the requirement in the interpolation theory thatNI(xJ)is equal to 1 whenI=Jand is equal to 0 whenI/=J.

Figure 9 : The convergence curves of the LS-based interpolation method with respect to the number of nodes Nn in example 1,the mean of relative variance error ∈Var (a),the mean of relative covariance error ∈Cov(b)

Figure 10 :The convergence curves of the LS-based interpolation method with respect to the number of nodes Nn in example 2,the mean of relative variance error ∈Var (a),the mean of relative covariance error ∈Cov(b)

Figure 11 :The convergence curves of the LS-based interpolation method with respect to the number of nodes Nn in example 3,the mean of relative variance error ∈Var (a),the mean of relative covariance error ∈Cov(b)

Table 1 displays other types of covariance functions with their corresponding parameters and domains.Meanwhile,the changes in the mean of the relative covariance error∈Covconcerning the number of truncating termsMare illustrated in Fig.12.The results indicate that the MLS-,FEMbased interpolation method and the FEM-based Galerkin method yield similar outcomes regarding∈Cov,with the error decreasing asMincreases.Conversely,the LS-based interpolation method produces comparable results to the previous three methods for covariance functions(b)and(d),but larger values for functions(a)and(c).Notably,the Legendre-based Galerkin method yields a higher value of∈Cov,which implies lower accuracy when compared to the other methods.

Table 1 : The covariance functions,σ=1

Figure 12 :Variation of the mean of relative covariance error ∈Cov with the number M of truncation.(a)Triangular,(b)brownian-bridge,(c)uniformly modulated,(d)modified exponential

Table 2 presents a comparison of the computation times for solving Fredholm integral eigenvalue problems using various covariance functions.The methods employ the same number of nodes,while the number of truncating terms is fixed atM=20.The computations were performed on a computer equipped with an i7-4790 CPU@3.60 GHz.The findings demonstrate that the LS-based interpolation method is the most efficient,followed by the MLS-and FEM-based interpolation methods,with the FEM-and Legendre-based Galerkin method ranking last in terms of computational time.Notably,the Galerkin method utilizing Legendre polynomials requires an order of magnitude more time than the other methods.Hence,the interpolation method offers better computational efficiency than the Galerkin method.The reason for this is that in the process of solving eigenvalues and eigenfunctions,elementClkof matrixCin the Galerkin method as shown in Eq.(18)is a double integral,while elementLI(xJ)of matrixLin the interpolation method as shown in Eq.(26)is a single integral.

Table 2 : The time costs(s)for solving Fredholm integral eigenvalue problems

4.3 Discretization of Two-Dimensional Random Fields and Stochastic Analysis

The following example is used to analyze the discretization of random fields in two-dimensional regular and irregular domains and is applied to stochastic analysis [37].According to the onedimensional analysis,we find that the interpolation method based on MLS and FEM can ensure both accuracy and relative stability.Therefore,in the following examples,we utilize MLS-and FEM-based interpolation methods for the discretization of random fields.

The provided Fig.13 depicts a rectangular plate with a height ofH=1.5 m and a width of=1 m.The plate is fixed at the bottom and subjected to a uniform tensile stress ofσ=1 Pa along its upper edge.This is a plane stress problem,and Poisson’s ratio isν=0.3.The modulus of elasticityEis supposed to be a random field with a mean value ofμE=10 Pa and a standard deviation ofσE=2 Pa.The covariance function is presented as follows:

The correlation length along thexdirection isbx=0.5 m,while along theydirection,it isby=0.75 m.The covariance function takes on a two-dimensional form of the exponential kernel function,which exhibits analytical eigenvalues and eigenfunctions[37].In KL expansion,the number of truncating terms isM=3.The Hermite polynomials are selected for the polynomial chaos expansion,and the orderp=3.Fig.14 displays the mesh used in the stochastic finite element method.

Figure 13 :A rectangular plate

Fig.15 illustrates the probability density function of the vertical displacementUyAat pointAon the plate.As seen that the probability density function calculated by the stochastic finite element method agrees well with the result by Monte Carlo(MC).Furthermore,the computed results based on the MLS and FEM techniques are found to closely match the analytical eigen-solutions,thus providing evidence for the validity of the MLS-and FEM-based interpolation methods for the KL expansion in two-dimensional domains.

Figure 15 :The probability density function fUyA(ξ)of the vertical displacement at point A

Fig.16 presents the standard deviationσUyof the vertical displacement of the plate obtained using various methods.It is evident from (a),(b),and (c) that the variation trend ofσUyderived through the MLS-and FEM-based interpolation techniques closely matches that of the analytical solution.Figs.16d and 16e illustrate the absolute errors in displacement standard deviations for the MLSand FEM-based interpolation methods relative to the analytical solution.The numerical outcomes reveal that both methods exhibit small absolute errors.However,the error associated with the MLSbased interpolation method is one order of magnitude lower than that of the FEM-interpolation method,indicating that the MLS-based interpolation approach is more accurate than the FEM-based interpolation method in this example.

Figure 16 : The standard deviation σUy of vertical displacement of the plate: (a) MLS,(b) FEMinterpolation,(c) analytical,(d) absolute error of MLS and analytical,(e) absolute error of FEMinterpolation and analytical

To investigate the effect of the standard deviation of elastic modulus and polynomial chaos order on the standard deviation of the vertical displacement of a plate,observation pointAwas selected.Fig.17 depicts the standard deviationσUyAof the vertical displacement of pointAfor polynomial chaos orderp=1,2,and 3 as (a),(b),and (c),respectively.The results demonstrate that as the standard deviationσEof elastic modulus E increases,σUyAgradually increases,with negligible differences among the results obtained using MLS-,FEM-based interpolation,and analytical eigen-solutions.Furthermore,as the standard deviation of the input random parameters increases,higher polynomial orders are required to achieve better accuracy.

To compare the applicability of different methods for the discretization of random fields in twodimensional irregular domains,a square plate with a pentagonal hole,as shown in Fig.18,was calculated and subjected to random analysis.The length of the plate isL=2 m,the left edge of the plate is fixed,while a uniformly distributed loadσ=1 Pa is applied to the other edge.The central coordinate of the hole is equal toxc=1,yc=1,the radius isR=0.3 + 0.08 × sin(5α),whereα∈[0,2π]rotates counterclockwise from the horizontal direction.The model is a plane stress problem,with Poisson’s ratioν=0.3.It is assumed that the elastic modulusEis a Gaussian random field,with a meanμE=20 Pa,and a standard deviationσE=5 Pa.The covariance function of the random field is the same as that of the previous example,and the correlation lengths arebx=1 m andby=1 m.

As the plate with a hole has an irregular domain,the basis function in the extended finite element method can be used to replace the finite element basis function when IEVP is solved in the discretization of random fields.The level set function is used to determine whether the node is in the plate.If the elements are cut by hole edge,only the part of the elements in the plate is integrated.The method is expressed as XFEM-interpolation.Meanwhile,in the subsequent stochastic analysis,the extended finite element method is combined with PCE,and the number of truncating terms isM=3,the order of polynomial chaos isp=3.

Figure 18 :A square plate and the XFEM mesh

An illustration of the different distribution modes of the MLS-based interpolation method used to solve the IEVP is presented in Fig.19.Specifically,Fig.19a shows a uniform distribution without considering the internal hole boundary,while Fig.19b shows a scattering distribution that takes into account the internal hole boundary.A comparison between the uniform and scattering distributions are conducted by the eigenvalues which are presented in Table 3.the results show that the errors of the two distribution mode are approximately 10-4.These findings suggest that different distribution modes have little impact on the accuracy of the discretization of random fields,and for convenience,the distribution mode shown in Fig.19a is selected for subsequent computations.

Figure 19 :The distributions of points of MLS-based interpolation method,(a)uniform distribution,(b)scattering distribution

Table 4 displays the standard deviationsσUxAandσUxBof the horizontal displacement at pointAand pointB,respectively,which were obtained by the MLS-based interpolation method under various correlation lengths.The table demonstrates that bothσUxAandσUxBgradually increase with increasing correlation lengths in both thexandydirections.Moreover,if one correlation length is fixed and the other increases,bothσUxAandσUxBalso gradually increase.Additionally,the impact of the correlation length in thexdirection on the standard deviation of the horizontal displacement is significantly greater than that of theydirection.

Table 4 :The standard deviation of horizontal displacement at points A and B with different correlation lengths

Fig.20 illustrates the distribution of the standard deviation of horizontal displacementσUxof the plate obtained using MLS-based interpolation method (a) and XFEM-based interpolation method(b).The discretization of random fields used 10×10 nodes based on MLS and 21×21 mesh nodes based on XFEM basis functions,while the stochastic analysis employed a mesh node of 21×21 for XFEM.The figure shows that the two methods produced similar trends in the distribution ofσUx,with the maximum value located at the middle of the right end of the plate.This suggests that the basis functions of XFEM can serve as interpolation basis functions for computing in complex domains.

Figure 20 :The standard deviation of horizontal displacement σUx:(a)MLS,(b)XFEM-interpolation

The XFEM calculations presented in the previous discussion employed the same mesh for both the discretization of the random field and stochastic analysis.However,the XFEM permits different meshes for the two processes.Nevertheless,applying eigenfunctions during stochastic analysis involves multiple transformations between global and local coordinates.In contrast,the MLS-based interpolation method enables the separation of distribution points and meshes in random field discretization and stochastic analysis without requiring any coordinate transformation.Hence,when a large number of meshes are used for stochastic analysis,the MLS-based interpolation method by point interpolation without meshing confers more advantages.

5 Conclusions

This study proposes an interpolation method for solving Fredholm integral eigenvalue problems in KL expansion for random field discretization.The performance of three interpolation basis functions,namely MLS,LS,and FEM,is evaluated.Compared to the Galerkin method,which uses finite element or Legendre polynomials as basis functions and involves a two-folded integral to calculate the integral matrix containing the covariance function,the proposed interpolation method only requires a single integral,resulting in reduced computational time.Numerical examples in one-dimensional domains reveal the validity and computational efficiency of the proposed method.The LS-based interpolation method is the most efficient,and the closer the multiplier of the node support domain is to 1,the higher the accuracy is.While the MLS-based interpolation method produces more stable results than LS,and the influence of the multiplier decreases with a large number of nodes.The FEM-based interpolation method has high accuracy,and the accuracy of the eigenfunction is almost independent of the eigenfunction index.

Additionally,this study combines KL expansion and PCE to perform random analysis in twodimensional regular and irregular domains.The results show that MLS-based interpolation provides higher computational accuracy than FEM-based interpolation in two-dimensional regular domains.As the standard deviations of the input random parameters increase,higher orders of polynomials are needed to achieve better results.In two-dimensional irregular domains,XFEM is used for stochastic analysis.The XFEM basis function can serve as the interpolation basis function,and MLS-based interpolation can be performed using uniformly distributed points for random field discretization.

Acknowledgement:The authors are grateful for the support by the Postgraduate Research&Practice Program of Jiangsu Province(Grant No.KYCX18_0526),the Fundamental Research Funds for the Central Universities (Grant No.2018B682X14) and Guangdong Basic and Applied Basic Research Foundation(No.2021A1515110807).

Funding Statement: The authors gratefully acknowledge the support provided by the Postgraduate Research & Practice Program of Jiangsu Province (Grant No.KYCX18_0526),the Fundamental Research Funds for the Central Universities (Grant No.2018B682X14) and Guangdong Basic and Applied Basic Research Foundation(No.2021A1515110807).

Author Contributions:The authors confirm contribution to the paper as follows:study conception and design:Zi Han;data collection:Zi Han;analysis and interpretation of results:Zi Han;draft manuscript preparation:Zi Han,Zhentian Huang.All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials:The data undering this article will be shared on reasonable request to the corresponding author.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.