Semi-reduced order stochastic finite element methods for solving contact problems with uncertainties

This paper develops two-step methods for solving contact problems with uncertainties. In the first step, we propose stochastic Lagrangian multiplier/penalty methods to compute a set of reduced basis. In the stochastic Lagrangian multiplier method, the stochastic solution is represented as a sum of products of a set of random variables and deterministic vectors. In the stochastic penalty method, the problem is divided into the solutions of non-contact and possible contact nodes, which are represented as sums of the products of two different sets of random variables and deterministic vectors, respectively. The original problems are then transformed into deterministic finite element equations and one-dimensional (corresponding to stochastic Lagrangian multiplier method)/two-dimensional (corresponding to stochastic penalty method) stochastic algebraic equations. The deterministic finite element equations are solved by existing numerical techniques, and the one-/two-dimensional stochastic algebraic equations are solved by a sampling method. Since the computational cost for solving stochastic algebraic equations does not increase dramatically as the stochastic dimension increases, the proposed methods avoid the curse of dimensionality in high-dimensional problems. Based on the reduced basis, we propose semi-reduced order Lagrangian multiplier/penalty equations with two components in the second step. One component is a reduced order equation obtained by smooth solutions of the reduced basis and the other is the full order equation for the nonsmooth solutions. A significant amount of computational cost is saved since the sizes of the semi-reduced order equations are usually small. Numerical examples of up to 100 dimensions demonstrate the good performance of the proposed methods.


Introduction
Contact problems play important roles in many practical applications [26,27].However, uncertainties are unavoidable for many problems in practical engineering.The consideration of uncertainties in systems behavior has become an essential part of the analysis and design of engineering applications.Although mathematical theories and numerical methods for the treatment of contact problems are well understood for the analysis of deterministic problems [6,9,27,28], there is only limited experience in the treatment of stochastic contact problems.
This paper focuses on solving two-body contact problems with random parameters, e.g.random material properties.There are several methods for solving stochastic contact problems.The Monte Carlo method (MC) is a powerful tool for many kinds of stochastic problems.Although MC can be readily implemented by use of the already existing deterministic codes, a large number of deterministic contact problems need to be solved in order to achieve high-accuracy solutions, which is prohibitively expensive for large-scale problems in practice [7].As an extension of MC, the multilevel Monte Carlo method (MLMC) has been developed to save the computational cost of MC.It has been applied to the model of contact between a deformable body and a rough uncertain substrate in [5] and the effect of uncertainties on quantities of interest related to contact mechanics of rough surfaces in [24].The adaptive MLMC is also presented to solve the stochastic obstacle problem in [16].
Another approach is the polynomial chaos (PC)-based method [12,29].In this method, the stochastic solution is projected onto a stochastic space spanned by (generalized) PC basis.The stochastic Galerkin projection is used to transform the original stochastic problem into a system of deterministic equations whose size is much larger than the original stochastic problem.However, the accuracy and efficiency of this method need to be further improved.On one hand, since the computational effort increases dramatically as the number of random variables and/or the number of expansion terms increase, the PC-based method suffers from the curse of dimensionality [2], which is prohibitively expensive, especially for large-scale and high-dimensional stochastic problems.On the other hand, the classical PC-based method cannot accurately capture the discontinuous solutions of stochastic contact problems, thus it usually works in conjunction with other methods, e.g. the stochastic collocation method [2,3,10].Other methods have also been proposed to solve related problems.The stochastic variational inequality problem (SVIP) is solved by stochastic projections and reformulations of SVIP in [14].Low-rank tensor methods are used to solve parametrized SVIP of the obstacle type in [11].
As described so far, efficient methods for solving stochastic contact problems are still insufficient and further investigations are needed.In this paper we propose efficient two-step algorithms to solve (high-dimensional) stochastic contact problems.In the first step, the stochastic Lagrangian multiplier method and the stochastic penalty method are developed to compute a set of reduced basis.The stochastic solution in the stochastic Lagrangian multiplier method is expanded as a sum of products of a set of random variables and deterministic vectors.The stochastic solution in the stochastic penalty method is divided into the solutions of non-contact and possible contact nodes, which are represented as sums of the products of two different sets of random variables and deterministic vectors, respectively.Similar expansions of stochastic solutions have been widely used in the proper generalized decomposition (PGD)/generalized spectral decomposition (GSD) methods [19,22].In this kind of method, priori tensor product constructions and separated variables representations are adopted to approximate the solutions.The components of the approximations are calculated by several dedicated iterations [8,17].PGD-based methods have been widely applied to various problems with great success [8,25] and the applications on contact problems can be found in [13,18].
The improvements of the approximations used in this paper include two aspects: one is that the non-intrusive representation of expanded random variables allows complexly nonlinear stochastic problems to be dealt with in a weak-intrusive way and to avoid the curse of dimensionality successfully, and the other is that a dedicated expansion is developed for the stochastic penalty method, which overcomes the poor convergence of solutions using the same extended random variables for both non-contact and possible contact nodes.
Based on the two proposed expansions, original stochastic finite element equations are transformed into deterministic finite element equations and one-dimensional (corresponding to the stochastic Lagrangian multiplier method)/twodimensional (corresponding to the stochastic penalty method) stochastic algebraic equations including all random variables.They are computed in their individual spaces, that is, deterministic finite element equations are solved by existing finite element methods and corresponding one-/twodimensional stochastic algebraic equations are solved by a sampling method, the computational cost of which does not increase dramatically as the stochastic dimension increases.Hence the curse of dimensionality in high-dimensional stochastic spaces is circumvented with great success.The effectiveness of the proposed approaches is demonstrated on numerical examples of up to 100 dimensions.
As the so far obtained full-reduced order equations directly constructed by the reduced basis cannot capture the nonsmooth solutions of stochastic contact problems immediately, semi-reduced order Lagrangian multiplier/penalty equations are proposed to overcome that obstacle.This strategy consists of the reduced order model into two components.One component is a reduced basis model obtained by smooth solutions and the other is the full order equation for the nonsmooth solutions.The size of the semi-reduced order equation is small since the possible contact region is small compared to the whole domain in many applications, thus the computational effort is significantly saved.Similar ideas have been discussed for the static condensation process of deterministic contact problems [21,26] and problems with localized solution behavior [1,15], in which the reduced approximation basis is adopted in the region with smooth solutions and a standard finite element approximation is used in the region with high gradient or discontinuous solutions.In this paper, the reduced basis in the first step is still considered as global approximations and smooth parts of the reduced basis are used to generate the semi-reduced order equations in the second step.
The paper is organized as follows: The basic setting of stochastic contact problems is introduced in Sect. 2. In Sect.3, we propose the stochastic Lagrangian multiplier method to solve stochastic contact problems, including the construction of stochastic solutions, the iterative algorithm and the algorithm implementation.Following that, a stochastic penalty method is developed in Sect. 4. In Sect.5, two numerical examples of low-and high-dimensional cases are given to demonstrate the performance of the proposed methods.Section 6 gives conclusions and discussions.

Stochastic contact problems
In this paper, we consider linear elastic two-body normal contact problems without friction and adopt the finite element method to discretize the problem.The total energy of the deterministic system is given by where U , W are the strain and kinetic energy and G comes from the contribution of the contact traction.Let ( , , P) be a complete probability space, where denotes the space of elementary events, is the σ -algebra defined on and P is the probability measure.Considering random parameters θ ∈ and extending Eq. ( 1) to the stochastic case, the total energy of the stochastic system is described as where U (θ ), W (θ ) and G (θ ) are considered as the stochastic strain energy, the stochastic kinetic energy and the stochastic contact contribution.They have similar forms to the deterministic cases but are related to the random input θ .The solution of Eq. ( 2) is obtained by minimizing the total potential energy (θ ), that is Several methods can be used to formulate G (θ ), e.g.Lagrangian multiplier method, penalty method and augmented Lagrangian method, etc [26].In this paper we only focus on the Lagrangian multiplier method and the penalty method.Other methods are beyond the scope of this article and will be studied in subsequent research.Similar to the deterministic case, G L (θ ) of the stochastic Lagrangian multiplier method (SLMM) is given by where s (θ ) is the stochastic contact interface, γ (θ ) is the vector of stochastic Lagrangian multipliers and coincides with the normal contact traction, g (θ ) is the random gap vector.
Corresponding the first variation δG L (θ ) in Eq. ( 3) is Similar to the classical penalty method, the stochastic G P (θ ) of stochastic penalty method (SPM) is given by where α is the penalty parameter and the corresponding first variation δG P (θ ) in Eq. ( 3) is calculated as Analogous to the deterministic cases [26], the advantage of SLMM is that it meets contact constraint conditions exactly.
The drawback is that finite element matrices based on SLMM are not positive definite and have augmented sizes due to introducing additional variables.It may present difficulties in dealing with some cases with friction.SPM has no additional variables but leads to very small physical penetration.We will develop numerical algorithms for solving the stochastic solutions based on SLMM and SPM in subsequent sections.

Stochastic Lagrangian multiplier method
In this section, we solve stochastic contact problems by adopting SLMM.By use of Eqs.(3) and ( 4) and the classical finite element discretization, a stochastic finite element equation regarding SLMM is obtained where K (θ ) ∈ R n×n is the stochastic stiffness matrix, F(θ ) ∈ R n is the stochastic force vector and Q ∈ R n c ×n is the matrix related to n c contact constraints and Lagrangian multipliers.The details of above finite element discretization can be found in [26].
It is noted that the coefficient matrix of Eq. ( 8) is not positive definite since is indefinite.In order to avoid this difficulty, we rewrite Eq. ( 8) in the form where The stochastic matrix K (θ ) in Eq. ( 11) is positive definite since x T , y T K (θ ) holds.In the next subsection, to take advantage of the positive definiteness of the matrix K (θ ), we solve the reduced basis based on Eq. ( 10) rather than Eq. ( 8).

Stochastic solutions of SLMM
To efficiently solve the stochastic solution u (θ ) of Eq. ( 10), we introduce the following approximation where the deterministic vectors )×k , the random variable vector (θ ) = [λ 1 (θ ), . . ., λ k (θ )] T ∈ R k and the previous approximation u k−1 (θ ) is given by For the stochastic components u(θ ) and γ (θ) of u (θ ), Eq. ( 13) is equivalent to where the vectors d i ∈ R n , q i ∈ R n c and each pair (d i , q i ) shares the same random variable coefficient λ i (θ ).It is noted that all of them are unknown and iterative algorithms are used to compute them one by one.The expansion Eq. ( 13) has been widely used in the context of PGD/GSD methods [8,19,22,32].We develop improved expansions in this paper to deal with stochastic contact problems and high-dimensional random variables involved.
Based on Eq. ( 13), Eq. ( 10) is written as where u k−1 (θ ) is assumed to be known and the couple {λ k (θ ) , d k } is the unknown solution to be solved.However, it is not easy to determine the random variable λ k (θ ) and the vector d k simultaneously.We adopt the following iterative algorithm [8,22] to solve λ k (θ ) and d k : If the random variable λ k (θ ) has been known (or given an initial value), the vector d k is solved via the following deterministic finite element equation obtained by stochastic Galerkin method [12,29], where E {•} is the expectation operator and which can be rewritten as a compact form, where the deterministic matrix and vector are given by Equation ( 17) can be considered as a finite element equation for deterministic contact problems, thus existing contact FEM solver can be used to solve it efficiently [26,27,34].Recalling Eq. ( 13), we approximate the stochastic solution by random linear combinations of the deterministic basis functions d k .Repeated or close modes of d k make the number of retained terms for the stochastic solution large.Therefore, in practice, we only want to retain the vector d k in dominant modes, so that the number of retained terms is as small as possible, which will also reduce the size of subsequent reduced order stochastic finite element equations.For this purpose, we let the solution d k orthogonal to the obtained vectors If the deterministic vector d k has been solved by use of Eqs. ( 17) and ( 20), the random variable λ k (θ ) is updated via the following deterministic Galerkin procedure where the simplification (12).A commonly used method to solve Eq. ( 21) is the polynomial chaos (PC)-based method, which expands λ k (θ ) by use of PC basis and solves the unknown expanded coefficients using a system of linear equations.However, this kind of method suffers from the curse of dimensionality when dealing with high-dimensional problems.It is noted that Eq. ( 21) is a one-dimensional stochastic algebraic equation of the random variable λ k (θ ).In this paper, we propose a sampling method to compute λ k (θ ) [32,33], which is corresponding to, where θ ∈ R × ×n s is the sample representation of stochastic matrix (or vector, variable) (θ ) ∈ R × , n s is the number of random samples, the vector represents the random sample vector of λ k (θ ).In practical calculations, Eq. ( 22) is performed by the element-wise division of the two vectors, which is similar to the Hadamard division operator.Equation ( 22) is less sensitive to the stochastic dimension and the computational effort is very low even for high stochastic dimensions, which avoids the curse of dimensionality successfully.Further, statistical methods can be used to compute probability characteristics of λ k (θ ) from the random samples λ k θ .
In this way, we calculate the k-th couple λ k (θ ) , d k by iteratively solving Eqs. ( 17) and ( 22) until a specified precision is reached, the details of which will be discussed in Algorithm 1. Also, the same iteration is used to calculate the next couple λ k+1 (θ ) , d k+1 until all couples are obtained.However, it is noted that the above iteration solves each couple λ k (θ ) , d k in a sequential way and the stochastic solution u (θ ) in Eq. ( 13) is also approximated in a greedy way.The original stochastic finite element equation (10) is not exactly met by the stochastic solution u k (θ ) obtained in this way.Low-accuracy stochastic solutions are obtained for some problems, as illustrated in Example 5.1.To improve the accuracy of the stochastic solution, we consider D = d 1 , . . ., d k as a set of reduced basis and introduce reduced order equations to recalculate the random vector

Full-reduced order method based on SLMM
As discussed above, we first build a full-reduced order equation by applying the known reduced order matrix D obtained by the above iteration to Eq. ( 10) which is equivalent to where the random vector (θ ) ∈ R k is the unknown solution to be recalculated, K D (θ ) = D T K (θ ) D ∈ R k×k is the reduced order stochastic matrix and R k is the reduced order stochastic force vector.In practice, we solve Eq. ( 24) by for each random sample realization θ (i) , i = 1, . . ., n s and the random samples θ = θ (1) , . . ., θ (n s ) ∈ R k×n s are obtained by n s solutions of Eq. (25).It is noted that the random samples θ = θ (i) n s i=1 ∈ used for solving Eq. ( 25) are not necessary to be the same as that used for solving the reduced order matrix D since only the deterministic matrix D is inherited from the iteration and no random quantities are involved.
The size of Eq. ( 24) is usually small thanking to the small value k, thus the computational effort for solving a single Eq. ( 25) is very low and the total computational costs of n s solutions are still low.However, the full-reduced order equation (24) cannot accurately capture the nonsmooth stochastic solutions on the contact interface, which will be verified in Example 5.2.Hence, it is necessary to improve the accuracy of Eq. ( 24).

Semi-reduced order method based on SLMM
In order to improve the computational accuracy of the stochastic solution of the full-reduced order equation and to capture the nonsmooth stochastic solution on the stochastic contact interface, we propose a semi-reduced order method to solve stochastic contact problems.
Recall the full order stochastic finite element equation ( 8) and rewrite it in the form, where the vectors u 1 (θ ) ∈ R n−n c and u 2 (θ ) ∈ R n c are the stochastic solutions of non-contact and possible contact degrees of freedom, the deterministic matrices ) and B 2 ∈ R n×n c are used to locate the non-contact and possible contact solutions and their elements are given by where δ i j is the Kronecker delta that meets δ i j = 1 for i = j and 0 for i = j, the set c con = c con,1 , . . ., c con,n c ∈ R n c contains the indices of the degrees of freedom on the possible contact interface and the set c non = c non,1 , . . ., c non,n−n c ∈ R n−n c are the indices of the rest of degrees of freedom.Let the submatrix )×k ⊂ D be the reduced basis matrix corresponding to the noncontact nodes.The stochastic solution of the non-contact nodes is thus represented as T ∈ R k is the random vector to be solved.Hence, the solution vector of Eq. ( 26) can be rewritten as where I n c ∈ R n c ×n c is the identity matrix.Based on Eq. ( 26), we introduce the following deterministic Galerkin procedure to build a semi-reduced order stochastic finite element equation which is rewritten as a compact form where the reduced order matrices and vectors are given by It is noted that only the stochastic solution of non-contact degrees of freedom is approximated by a reduced order representation and the stochastic solution of possible contact degrees of freedom is still solved in a full order way, thus the proposed method is a semi-reduced order approach.Further, the size of the semi-reduced order equation ( 31) is k + 2n c , which is slightly larger than the size k of Eq. ( 24), but still very small compared to the original stochastic problem Eq. ( 8) in most applications.Similar to Eq. ( 25), we still adopt the sampling method to solve (31) for each sample realization and total computational costs of n s solutions are still very low.

Algorithm implementation of SLMM
The proposed SLMM for solving stochastic contact problems is summarized in Algorithm 1, which includes a doubleloop iterative procedure.The inner-loop iteration from step 3 to 10 is used to compute the couple λ k (θ ) , d k and the outer loop from step 2 to 15 corresponds to recursively approximate the stochastic solution u k (θ ).In order to compute each couple λ k (θ ) , d k , the random sample vector In practical implementation, any nonzero vectors with size n s can be used as the initialization and it almost has no influence on the proposed method.The vector d k, j is determined by solving the deterministic finite element equations in step 5, where the subscript j represents the j-th round iteration of λ k, j (θ ) , d k, j .Following that, d k, j requires to be orthogonalized and normalized along the whole iterative process.With the obtained vector d k, j , the random variable λ k, j (θ ) is then updated in step 7.The iterative error in step 8 is given by where = E T .The outer-loop iteration then generates a new approximated stochastic solution u k (θ ) in step 11 and updates the deterministic matrix D in step 12.The Algorithm 1 SLMM for solving stochastic contact problems 1: Discretize and generate the stochastic finite element equation ( 8) 2: while ε g,k > ε g do 3: Initialize random samples while ε l, j > ε l do 5: Compute d k, j by solving Eq. ( 17) i=1 by using Eq. ( 20) and normalize d k, j = 1 7: Compute λ k, j θ ∈ R ns via Eq.( 22) 8: Compute the iterative error ε l, j 9: j ← j + 1 10: end while 11: Update the stochastic solution Update the deterministic matrix D = D, d k ∈ R (n+nc)×k 13: Compute the iterative error ε g,k 14: k ← k + 1 15: end while 16: Solve the semi-reduced order stochastic finite element equation (31) iterative error in step 13 is defined as which measures the contribution of the k-th stochastic increment λ k (θ ) d k to the stochastic solution u k (θ ).Based on the known matrix D, the semi-reduced order stochastic finite element equation with the size k + 2n c is obtained and solved in step 16.We highlight that the sample size n s used for solving the reduced basis can be different from that for solving the semi-reduced order stochastic equation.In practice, we can adopt a small size in step 3 to calculate the reduced basis and after obtaining the matrix D, a large sample size is reset in step 16 to solve the semi-reduced order stochastic finite element equation, which saves computational costs.

Stochastic penalty method
In this section, we propose a SPM to solve stochastic contact problems.According to Eqs. ( 3) and ( 6) and the classical finite element discretization, the stochastic finite element equation obtained by the penalty method is described as where K 0 (θ ) ∈ R n×n is the stochastic stiffness matrix, α is the penalty parameter, K p ∈ R n×n is a deterministic matrix related to possible contact constraints and the detail of its assembly can be found in [26].
In the approximation Eq. ( 13), the vectors d k of possible contact and non-contact degrees of freedom share the same random variable coefficients {λ k (θ )}.It does not work well for SPM since a large penalty parameter biases the random variables {λ k (θ )} towards the solution of the possible contact degrees of freedom.Let us give a simple explanation for this point.For the purpose, we still consider an Eq. ( 13)-like approximation used to SPM, that is, the coefficient of the deterministic vector which is further rewritten as the approximate equivalence in which holds for a large penalty parameter α.The penalty parameter biases λ * k (θ ) towards the contribution of d * T k K p d * k , which cannot capture the contribution of non-contact nodes well.

Stochastic solutions of SPM
In order to solve Eq. ( 35), we introduce the following approximation to avoid the difficulty suffering from the expansion Eq. ( 13), where the vector u k−1 (θ ) ∈ R n is the previous approximation of the stochastic solution and it is usually already known, the vector u k (θ ) ∈ R n is the unknown stochastic increment to be solved.The size of u k (θ ) in Eq. ( 38) is the same as the original stochastic problem, which is different from the augmented vector u k (θ ) ∈ R n+n c in SLMM.
To calculate the unknown stochastic increment u k (θ ), we introduce the following two equivalent approximations where d k ∈ R n is the unknown deterministic vector to be solved, the random variable vector 27) and (28).They can also be calculated by where the matrix Similar to SLMM, an alternating iteration is adopted in this section to solve the triplet λ k,1 (θ ) , λ k,2 (θ ) , d k in Eq. (39) one by one.Specially, for known random variables λ k,1 (θ ) and λ k,2 (θ ), the vector d k is solved via the following deterministic finite element equation obtained by the approximation Eq. (39) and stochastic Galerkin method 123 where the deterministic matrix and vector are given by Further, Eq. ( 42) can be rewritten as where the matrices K 0,k and K p,k are obtained by Eq. ( 43) In this way, Eq. (42) (or Eq. ( 45)) can also be considered as a finite element equation with penalty parameters for deterministic contact problems and existing contact FEM solvers are available.Similar to SLMM, we retain the dominant vector d k and Eq. ( 20) is still adopted to orthogonalize the vector d k .
Based on the obtained vector d k , an updated random vector T is solved by Eq. ( 40) and the following deterministic Galerkin procedure which is rewritten as where the random variables c k,i j (θ ) and h k,i (θ ) are given by for i, j = 1, 2. By taking advantage of the proposed sampling method used in Eq. ( 22), the stochastic solutions λ k,1 (θ ) and λ k,2 (θ ) of Eq. ( 49) are solved via the following two decoupled equations where the random sample vectors of the numerators s k,1 θ , s k,2 θ ∈ R n s and the denominator s k,0 θ ∈ R n s are given by Computational costs for solving Eqs. ( 52) and ( 53) are still very low.Also, similar to Eq. ( 22) used in SLMM, both Eq. ( 52) and ( 53) are insensitive to the stochastic dimension, thus SPM can also be applied to high-dimensional stochastic contact problems.

Semi-reduced order method based on SPM
In order to build a reduced order equation for SPM, we rewrite the approximate stochastic solution u k (θ ) as where is a random vector and u 2 (θ ) ∈ R n c is the stochastic solution of the possible contact degrees of freedom.Similar to Eq. ( 30), the semi-reduced order stochastic finite element equation of SPM is then given as where the sizes of the stochastic submatrices and subvectors are . The sampling approach used in Eq. ( 25) is again adopted to solve the semi-reduced order stochastic finite element equation (56) for each sample realization.One-sample solution of Eq. ( 56) is cheaply solved since its size k + n c is much smaller than the original stochastic problem.Total computational costs for n s random samples θ = θ (1) , . . ., θ (n s ) ∈ R (k+n c )×n s of the solution are still low.

implementation of SPM
Algorithm 2 SPM for solving stochastic contact problems 1: Discretize and generate the stochastic finite element equation (35) 2: while ε g,k > ε g do 3: Initialize the random samples λ (0) k by solving Eq. (42) 6: Orthogonalize d k, j ⊥ {d i } k−1 i=1 by using Eq. ( 20) and normalize d k,2 θ ∈ R ns via Eqs.( 52) and (53) 8: Compute the iterative error ε l, j 9: j ← j + 1 10: end while 11: Update the stochastic solution Compute the iterative error ε g,k 14: k ← k + 1 15: end while 16: Solve the semi-reduced order stochastic finite element equation ( 56) The proposed stochastic penalty method for solving stochastic contact problems is summarized in Algorithm 2 and the procedure is similar to Algorithm 1.It is noted that samples of two random variables λ k,1 (θ ) and λ k,2 (θ ) need to be initialized in step 3 and corresponding updates of two random variables are processed in step 7.After the innerloop iteration, the stochastic increment u k (θ ) in step 11 is updated via Eqs.(39) or (40).Similar to Eq. ( 34), the iterative error in step 13 is now defined by use of random variables and the iterative error ε l, j in step 8 has been defined in Eq. ( 33).

Numerical examples
In this section, we test the proposed SLMM and SPM by two numerical examples.For both SLMM and SPM, the convergence errors are set as ε l = 1 × 10 −3 and ε g = 1 × 10 −6 in Algorithm 1 and 2. The number of random samples is n s = 1 × 10 4 and the reference solutions are obtained by 1 × 10 4 standard MCS.The penalty parameter is α = 10 6 N/m 2 in SPM.
Fig. 1 Euler Bernoulli beam supported by a rigid block with a gap

Euler Bernoulli beam limited by a rigid block
This example considers a one-dimensional Euler Bernoulli beam shown in Fig. 1, which is clamped at one end and subjected to a distributed load q(x).The deflection of the free end of the beam is limited by a rigid block.There is a small random gap between the end of the beam and the rigid block.The bending rigidity E I (x, θ) is considered as a Gaussian random field with the mean function E I (x) and the covariance function given by By use of Karhunen-Loève (KL) expansion [4,31], the random field E I (x, θ) is approximated by the following series expansion where r is the truncation number, {ξ i (θ )} r i=1 are independent standard Gaussian random variables and {κ i , f i (x)} r i are solved by the homogeneous Fredholm integral equation of the second kind, In this example the following numerical parameters are adopted: the distributed load q(x) = 1 × 10 3 kN/m, the length of the beam L = 1 m, the mean value of the bending rigidity E I (x) = 1 × 10 7 N • m 2 , the standard deviation σ E I = 0.2E I (x), the correlation length l x = 1 m, the truncation number r = 5 and the random gap is a uniform random variable on [1, 1.5] × 10 −2 m.It is noted that the random samples θ (i) such that min x E I x, θ (i) ≤ 1 × 10 −3 are dropped out to make sure that all realizations are positive, thus E I (x, θ) is considered as a truncated Gaussian random field in the numerical implementation The model is discretized into 100 two-node Euler Bernoulli beam elements [23].In this element, each node includes 2 degrees of freedom given by the vertical displacement and the rotation.Thus Figure 2 shows the comparison between the probability density functions (PDFs) of the stochastic solution of the vertical displacement of the free end without the rigid block and the random gap.The free end is in the contact state for the stochastic solution exceeding 0.015 m and is in the noncontact state for the stochastic solution less than 0.01 m.The stochastic contact occurs for the solution in the random gap interval.In other words, the occurrence of the contact state between the free end and the rigid block is random since the stochastic solution and the gap are independent in the interval [1, 1.5]×10 −2 .It is necessary to check the stochastic contact state in practical computations.
We make use of the proposed Algorithm 1 and 2 to solve this problem.Iterative errors of the proposed SLMM and SPM are shown Fig. 3, which indicates that SPM is more efficient than SLMM in this case since only k = 3 retained items are needed to achieve the specified accuracy while SLMM needs nine items.
The PDFs of the stochastic vertical displacement of the free end obtained by SLMM, SPM and MCS are compared in Fig. 4, where the SLMM (initial) and SPM (initial) represent the solutions obtained by the iterative step 2 to 15 of Algorithm 1 and 2, i.e. the stochastic solution u k (θ ) in step 11 after the last iteration, and the SLMM (final) and SPM (final) represent the solutions obtained by semi-reduced stochastic finite element equations in step 16 of Algorithm 1 and 2. As shown in Fig. 4a, b, the solutions of both SLMM (final) and SPM (final) have good agreement with MCS, which demonstrates the high accuracy of the proposed methods.SLMM (initial) shows less accuracy and the semi-reduced stochastic finite element equation is necessary to improve the approximation accuracy while SPM (initial) is already of high accuracy.It is noted that the PDFs in Fig. 4 are discontinuous at the boundary values 0.01 m and 0.015 m of the random gap.In fact, as shown in Fig. 5, for each sample realization of the random gap, the PDF of the vertical displacement of the free end is discontinuous at the value of the sample.If all samples of the gap are considered in the case of Fig. 5, we can obtain the PDF that is continuous in the interval of the random gap but still discontinuous at the boundary.The proposed methods can capture this characteristic well, as depicted in Fig. 4.Moreover, PDFs of contact forces of the free end obtained by SLMM and MCS are shown in Fig. 6, which indicates the proposed SLMM is in good accordance with MCS and accurately captures the discontinuous contact force around zero (i.e. the non-contact state).

Contact between hemisphere and rigid body
In this example we consider a hemisphere in contact with a rigid body shown in Fig. 7a where the radius is 1.0 m and the force f (x, y) = 1.5 × 10 3 MN/m.We only use a half model for the analysis by taking advantage of the symmetry.The finite element mesh is depicted in Fig. 77b, including 374 nodes and 672 triangular elements with linear shape functions.Poisson's ratio of the hemisphere is ν = 0.2.The Young's modulus E (x, y, θ) is considered as a Gaussian random field with the mean value E(x, y) = 210 GPa and the covariance function where the standard deviation σ E = 0.1E(x, y) and the correlation lengths l x = l y = 1 m.Similar to Eq. (59), we approximate the random filed E (x, y, θ) by use of KL expansion.The number of the truncated item of KL expansion is adopted as r = 10.Again, the sample realizations θ (i) such that min x,y E x, y, θ (i) ≤ 1 × 10 −3 are dropped out to generate positive random samples.The iterative errors of SLMM and SPM are shown in Fig. 8.The items k = 11 for SLMM and k = 9 for SPM are retained to achieve the specified accuracy, which indicates that SPM is still more efficient than SLMM.In this sense, SPM is suggested to be a better choice for solving most stochastic contact problems.
The first six reduced basis in x and y directions obtained by SLMM and SPM are depicted in Fig. 9.As shown in Fig. 9a, b, the comparison of the reduced basis in the x direction between SLMM and SPM indicates that the reduced basis of possible contact nodes are more dominant components in SPM.It is seen from Fig. 9c, d that this conclusion is more obvious in the y direction since the contact boundary is closely related to the displacement in the y direction.Also, the first item d x1 and d y1 of SLMM are similar to that of SPM since both of them are solved by contactless stochastic problems.
To illustrate the accuracy of the proposed methods, we compare the solutions obtained by the proposed methods and MCS.As shown in Fig. 10, the solutions of five characteristic nodes are checked, including the contact nodes, the possible contact nodes and the non-contact nodes on the possible contact interface, the nodes near and far from the possible contact interface.For the contact nodes on the possible contact interface, SLMM matches contact boundary conditions accurately.The relative errors between the vertical displacements of contact nodes obtained by SPM and the gaps are shown in Table 1, where gap i represents the initial gap value of the i-th contact node (i.e. the cyan points) shown in Fig. 10.The results in Table 1 show the high accuracy of SPM for solving the solutions of contact nodes.The possible contact nodes (i.e. the red points in Fig. 10) on the possible contact interface are the nodes that are in contact states for some samples of random parameters and in non-contact states for other samples, thus they are in 'possible' contact states.Accurately capturing possible contact states is a difficulty for stochastic contact problems since their PDFs are strongly discontinu- We also check the nodes that are not on the possible contact interface, whose stochastic solutions are solved by the reduced order component of the semi-reduced order stochastic finite element equation.PDFs of the solutions of the points  A and B (as depicted in Fig. 10) in x and y directions are shown in Fig. 13.The good accuracy of SLMM and SPM is certified again.
For comparisons, we solve the stochastic solutions by use of the full-reduced order method proposed in section 3.2 and the PC-based stochastic Lagrangian multiplier method given in "Appendix".The size of the full-reduced order equation is equal to the retained items k = 11.Two-order Hermite PC basis of ten standard Gaussian random variables are adopted.The number of PC basis is 66 and the size of the derived finite element equations is about 5 × 10 4 , which is much higher than the original stochastic problem.Solutions of the contact nodes on the possible contact interface obtained by the full-reduced order method and PC-based SLMM are seen in Fig. 14, which indicates that both the full-reduced order method and PC-based SLMM cannot meet the contact boundary exactly.
where t S (second) is the computational time for computing the reduced basis, t U is the computational time for solving the semi-reduced order stochastic finite element equation and t T is the total time.
To verify the validity of the proposed methods for solving high-dimensional stochastic contact problems, different stochastic dimensions r = 10 and 100 are solved using SLMM and SPM.All of them are tested on a laptop (dualcore, Intel Core i7, 2.40GHz).Computational times are listed in Table 2, which indicates that the proposed SLMM and SPM are efficient even for the stochastic dimension 100.For both SLMM and SPM, as the stochastic dimension increases, the computational time t U for solving the semi-reduced order stochastic finite element equations slightly increases since the retained item k is a bit larger than the low-dimensional case while the computational time t S for computing reduced basis increases.In this paper, we adopt KL expansion (as shown in Eq. ( 59)) to approximate the random fields and the stochastic matrix K (θ ) in SLMM and SPM has the form K (θ ) = r i=0 ξ i (θ ) K i .Extra computational effort and storage are needed for a large number of deterministic matrices {K i } r i=0 , thus the computational time t S is high for the high-dimensional case.High-performance computational equipment and parallel computing can be used to reduce the computational time furthermore.Further, we study the influence of the sample size n s on the proposed methods and only SLMM used for the case of 10 stochastic dimensions is tested as an illustration.In numerical implementation, the sample sizes n s = 10 1 , 10 2 , 10 3 , 10 4 in step 3 of Algorithm 1 are first used to solve the reduced basis matrix D and n s is then reset as 10 4 in step 16 to solve the semi-reduced order stochastic finite element equation.PDFs of stochastic displacements of the point A in the y direction obtained by different sample sizes n s are compared in Fig. 15a and their absolute errors relative to the reference PDF obtained by MCS are seen from Fig. 15b, which indicates that the computational accuracy decreases slightly as the sample size n s decreases.It is noted that the PDF accuracy is acceptable in most cases even only 10 samples are adopted to calculate the reduced basis, thus we do not need to carefully choose the sample size in many problems.
Further, computational costs of different sample sizes are listed in Table 3, where the times t S , t U and t T have the same  meaning as that in Table 2.It is seen that for all cases, the computational times t U for solving the semi-reduced order stochastic finite element equations are almost the same since the number of retained terms is very close.While the computational times t S for calculating reduced basis decreases as the sample size decreases since the computational costs for performing Eqs. ( 18), ( 19) and ( 22) are reduced.

Conclusion
This paper develops two semi-reduced order stochastic finite element methods for solving stochastic contact problems and two numerical examples are used to certify their accuracy and efficiency.Based on the approximations of stochastic solutions proposed for SLMM and SPM and the corresponding iterative algorithms, the reduced basis are calculated by solving deterministic finite element equations and one-/two-dimensional stochastic algebraic equations.By use of the obtained reduced basis, the original problems are transformed into semi-reduced order stochastic finite element equations, whose solutions are solved by using a nonintrusive sampling method.The proposed methods can be applied to high-dimensional stochastic contact problems and avoid the curse of dimensionality successfully.In these senses, the proposed methods provide novel and efficient numerical strategies for solving stochastic contact problems.We only focus on the stochastic Lagrangian multiplier method and stochastic penalty method in this paper.Extending other methods for solving deterministic contact problems to the stochastic cases is attractive, e.g. the augmented Lagrangian multiplier method.Also, we only consider linear elastic contact problems.Contact problems coupled with other nonlinearities should be further studied, e.g.material nonlinearity or large deformations.
tive study between the proposed methods and the PC-based method in Example 5.2.In this section, we give a simple process to solve the stochastic finite element equation ( 8) using the PC-based stochastic Lagrangian multiplier method.
In the context of the PC-based method, the stochastic solution of Eq. ( 8) is represented as where { i (θ )} P i=1 are PC basis, the total number P of the PC basis is (r + p)!/r !/ p!, where (•)! represents the factorial operator, r and p are the number of random variables and the order of PC basis, respectively, the unknown vector u T i , γ T i T ∈ R n+n c are the deterministic expanded coefficients corresponding to the basis i (θ ).Although there are theoretical guidelines [29], the type of PC basis and the expansion order still need to be specified manually.Substituting Eq. (62) into Eq.( 8) and applying stochastic Galerkin method we have −g (θ ) j (θ ) , j = 1, . . ., P, (63) equivalently, where the matrices K ji = E K (θ ) i (θ ) j (θ ) ∈ R n×n , the vectors F j = E F (θ ) j (θ ) ∈ R n , g j = E g (θ ) j (θ ) ∈ R n c and which can be rewritten as (66) It is worth noting that the size of the PC-derived deterministic finite element equation ( 65) is n P = (n + n c ) (r + p)!/r !/ p!, which is significantly higher than the original stochastic problem.It is prohibitively expensive to solve Eq. (65) for high-dimensional and large-scale stochastic contact problems and thus leads to the curse of dimensionality.

Fig. 2 Fig. 3
Fig.2The solution of the free end without the rigid block and the random gap

Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8
Fig.4 PDFs of stochastic vertical displacement of the free end obtained by SLMM, SPM and MCS

Fig.Fig. 10 Fig. 11 Fig. 12
Fig. First six reduced basis {d xi } 6 i=1 in the x direction and d yi

Fig. 13
Fig.13 PDFs of the stochastic displacements of the points A and B

Fig. 14
Fig.14 PDFs of solutions of the contact nodes the possible contact interface obtained by full-reduced order method and PC-based SLMM

Fig. 15
Fig. 15 PDFs of stochastic displacements of the point A in the y direction obtained by and different sample sizes n s (top) and their absolute errors relative to the reference PDF obtained by MCS (bottom) ) where the solutions u = u T 1 , . . ., u T P T ∈ R n P , γ = γ T 1 , . . ., γ T P T ∈ R n c P and the deterministic matrices are given by

Table 3
Computational costs of different sample sizes n s