A nonlinear stochastic finite element method for solving elastoplastic problems with uncertainties

This article presents an efficient nonlinear stochastic finite element method to solve stochastic elastoplastic problems. Similar to deterministic elastoplastic problems, we describe history‐dependent stochastic elastoplastic behavior utilizing a series of (pseudo) time steps and go further to solve the corresponding stochastic solutions. For each time step, the original stochastic elastoplastic problem is considered as a time‐independent nonlinear stochastic problem with initial values given by stochastic displacements, stochastic strains, and internal variables of the previous time step. To solve the stochastic solution at each time step, the corresponding nonlinear stochastic problem is transformed into a set of linearized stochastic finite element equations by means of finite element discretization and a stochastic Newton linearization, while the stochastic solution at each time step is approximated by a sum of the products of random variables and deterministic vectors. Each couple of the random variable and the deterministic vector is also used to approximate the stochastic solution of the corresponding linearized stochastic finite element equation that can be solved via a weakly intrusive method. In this method, the deterministic vector is computed by solving deterministic linear finite element equations, and corresponding random variables are solved by a non‐intrusive method. Further, the proposed method avoids the curse of dimensionality successfully since its computational effort does not increase dramatically as the stochastic dimensionality increases. Four numerical cases are used to demonstrate the good performance of the proposed method.


INTRODUCTION
Since the inherent uncertainty and the epistemic uncertainty of systems are unavoidable in many practical engineering problems, 1,2 predicting uncertainty propagation on the physical models has become an important part of the analysis of systems. In this article, we focus on efficiently solving elastoplastic problems with uncertainties, in which the uncertainties are described by random variables or fields in the context of a probabilistic frame. Although mathematical theories and numerical methods for solving deterministic elastoplastic problems have been well understood, 3,4 only a few effort is took to their extensions to stochastic cases. Over the last two decades, several methods can be used for this purpose. As a kind of powerful method for stochastic analysis, Monte Carlo simulation (MCS) and its improvements 1,[5][6][7] can be used to accurately solve stochastic elastoplastic problems. This method is easy to be implemented and only requires the use of existing deterministic solvers. Its convergence does not depend on the stochastic dimensionality, thus MCS can be applied to solve high-dimensional stochastic problems without extra difficulties. However, a large number of deterministic simulations are required to achieve high-accuracy stochastic solutions, which is expensive for large-scale and nonlinear stochastic problems. Due to the simple and efficient implementation, other non-intrusive methods are also developed to solve stochastic elastoplastic problems, for example, adaptive sampling methods, surrogate model methods and so forth. [8][9][10] As a kind of intrusive method, the spectral stochastic finite element method (SFEM) and its extensions 1,11,12 have been widely used to solve stochastic problems. In this method, the stochastic solution is expanded by polynomial chaos (PC) bases and the stochastic Galerkin approach is then adopted to transform the original stochastic problem into an augmented deterministic equation whose size dramatically increases as the numbers of the stochastic dimensionality, the expansion order of PC bases and the degree of freedom of the physical system increase. Several improvements 13,14 are developed to reduce the computational effort of the augmented deterministic equation. A direct extension of this method to stochastic elastoplastic problems can be found in Reference 15. Several modifications are studied since the classical PC expansion cannot well capture nonsmooth stochastic solutions of stochastic elastoplastic problems. References 16 and 17 introduce two fictitious bounding bodies to approximate stochastic elastoplastic constitutive equations and then PC-based expansions are adopted to solve the corresponding stochastic finite element equations. Stochastic elastoplastic problems are solved by coupling the spectral SFEM to the Fokker Planck Kolmogorov (FPK) equation approach in References 18 and 19. The randomness of the elastoplastic constitutive model is propagated through FPK equations and the unknown stochastic solution is solved via the spectral SFEM. In Reference 20, stochastic elastoplastic problems are considered as inequality-constrained stochastic boundary value problems that are investigated via stochastic variational inequalities. PC expansions and stochastic collocation methods are used to discretize the original stochastic problem and nonlinear convex programming is adopted for the numerical implementation. Also, since the multi-element PC expansion 21 can be used to well capture sharp and discontinuous stochastic solutions, it has been applied to stochastic elastoplastic problems in Reference 22. There are also other effective methods for solving nonsmooth stochastic solutions of nonlinear stochastic problems, for example, the stochastic perturbation method, 23,24 the stochastic collocation methods, 25,26 the polynomial/spline dimensional decomposition methods. [27][28][29] As discussed above, developing efficient and accurate methods for solving stochastic elastoplastic problems remains an attractive topic. In this article, we propose a novel nonlinear SFEM to solve stochastic elastoplastic problems. Specifically, stochastic elastoplastic problems are considered as (pseudo) time-dependent stochastic nonlinear problems. A time discretization scheme is given to evolve the stochastic elastoplasticity and discretize the original time-dependent stochastic nonlinear problems into the time-independent stochastic nonlinear problem at each time step, and the finite element method is adopted for the spatial discretization. For each time step, the initial values of the discretized stochastic nonlinear problem are inherited from the stochastic quantities of the previous time step, for example, stochastic strains and stochastic internal variables. A stochastic Newton method is given to linearize the discretized stochastic nonlinear problem as a series of linear stochastic finite element equations. The stochastic solution at each time step thus consists of a set of stochastic increments that are also the stochastic solution of the corresponding linear stochastic finite element equation. Each stochastic increment is then decoupled into stochastic and deterministic spaces and approximated by the product of a random variable and a deterministic vector, which are calculated by solving the corresponding linear stochastic finite element equation using a weakly intrusive SFEM. 30,31 In this method, the deterministic vector is computed by solving deterministic linear finite element equations that are obtained by the stochastic Galerkin method, 11,12 and the corresponding random variable is solved by one-dimensional stochastic algebraic equations using a non-intrusive method. Furthermore, since all random sources are embedded into one-dimensional stochastic algebraic equations, the proposed method has low computational complexity and its computational cost is weakly dependent on the stochastic dimensionality. The curse of dimensionality arising in high-dimensional stochastic elastoplastic problems is thus circumvented successfully.
The article is organized as follows: Section 2 gives the basic setting of stochastic elastoplastic problems and the time discretization scheme. Section 3 introduces a stochastic Newton-based nonlinear stochastic finite element method to solve stochastic elastoplastic problems. Computational aspects and algorithm implementations of the proposed method are elaborated in Section 4. Following that, four numerical cases are given to demonstrate the performance of the proposed method in Section 5, and conclusions and discussions are presented in Section 6.

Elastoplastic equation with uncertainties
Let (Θ, Ξ, ) be a suitable probability space, where Θ denotes the space of elementary events, Ξ is a -algebra defined on Θ and  is a probability measure. In this article, we consider the following stochastic elastoplastic problem defined on a deterministic domain Ω ⊂ R d with the boundary Ω, where the physical dimension may be d = 1, 2, 3, div (⋅) represents the divergence operator, (x, ) is the stochastic stress tensor, f (x, ) is the stochastic external force, u (x, ) is the unknown stochastic displacement, Γ N and Γ D are boundary segments associated with the Neumann boundary condition g (x, ) and the Dirichlet boundary condition u (x, ). Without loss of generality, we let the Dirichlet boundary condition u (x, ) = 0 in this article. The following linear stochastic strain tensor is adopted to describe the stochastic displacement-strain relationship Further, we consider the relationship between the stochastic stress tensor and the stochastic strain tensor given by where  ∶ R d×d → R d×d is a stochastic nonlinear operator between the stochastic strain and the stochastic stress tensors. It may be associated with stochastic yield stresses and stochastic constitutive relations, such as the Von Mises and Drucker-Prager yield criteria 4 involving the stochastic thermodynamical force ( ) and the stochastic internal variable ( ). The stochastic operator  in this section is a general form and it can be applied to all kinds of elastoplastic models. We will provide detailed representations of it for some cases in the subsequent section.
The weak form of (1) is written as the left and right sides of which are given by where v (x) is the test function and  is the Hilbert space with smooth displacement vectors vanishing on the boundary Γ D .

Elastoplastic constitutive model
In this section, we consider applying an elastoplastic constitutive model to (1). To this end, we divide the stochastic strain tensor (u (x, )) into an elastic part and a plastic part, which corresponds to where e (u (x, )) and p (u (x, )) denote the stochastic elastic strain tensor and the stochastic plastic strain tensor, respectively. The stochastic stress tensor (u, ) is represented by the stochastic elastic strain tensor e (u (x, )) where C ( ) is a fourth order elastic tensor and the following form is adopted in this article where the parameters ( ) and ( ) are random variables/fields that may be related to stochastic material properties, also known as stochastic bulk modulus and stochastic shear modulus, the second order unit tensors I ∈ R d×d , the fourth order tensor I dev = I 4 − 1 3 I ⊗ I, where the fourth order tensor I 4 = ik jl e i ⊗ e j ⊗ e k ⊗ e l , ik is the Kronecker delta. Further, we only consider the Von Mises yield criterion and the linear kinematic hardening in this article, but other cases can also be considered in a similar way, for example, the Drucker-Prager yield criterion, the isotropic hardening and so forth. 4 The stochastic thermodynamical force ( ) associated with the internal variable ( ) is given by where the coefficient ( ) is a positive parameter and can be a random variable/field in practice. Furthermore, the stochastic yield function is given by where Y ( ) is the stochastic yield stress that can be modeled by a random variable/field, I dev ( ) is considered as the stress deviator of the stochastic stress tensor (u (x, ) , ).

Time discretization scheme
In this section, we extend the incremental method for solving deterministic elastoplastic problems to stochastic cases. The external force is applied incrementally in steps, which is considered as a set of pseudo time steps and described using a pseudo time variable t.
. For the time step t i , i ≥ 1, the stochastic solution u (t i , ) is calculated by solving a nonlinear stochastic problem that is dependent on the quantities } obtained by the previous time step t i−1 . For simplicity, we denote (t i , ), p (t i , ) and (t i , ) as i ( ), p,i ( ) and i ( ) in the latter equations. According to the analytical results in Reference 32, an explicit expression of the elastoplastic constitutive relation (3) can be obtained in the stochastic case, where i ( ) = 1 − Y ( ) | * i ( )| , the operator ⟨⋅⟩ = max(⋅, 0) represents the maximum value of ⋅ and 0, and the stochastic stresses 0,i ( ) and * i ( ) are given by The i ( ) and p,i ( ) are then updated by using * i ( ) from current time step t i and i−1 ( ), p,i−1 ( ) from the previous time step t i−1 , It is seen that ⟨ i ( )⟩ = i ( ) > 0 holds and the second terms of (13), (16), and (17) are nonzero when , the true stress is greater than the yield stress), the corresponding nodes are thus in plastic states. Otherwise, these second terms are all zero and we do not need to update (13), (16), and (17). Also, the stochastic tangent tensor  T of  is usually required to solve the linearized equation of the nonlinear problem at each time step. The analytical expression of  T is given by

STOCHASTIC NEWTON-BASED NONLINEAR STOCHASTIC FINITE ELEMENT METHOD
By use of the above relations and the finite element method for spatial discretization, we can obtain a random parameter-dependent nonlinear finite element equation, also known as the nonlinear stochastic finite element equation. In this section, we give an efficient and accurate numerical algorithm to solve nonlinear stochastic finite element equations.

Stochastic Newton linearization
For simplicity, we write the nonlinear stochastic finite element equation obtained by the elastoplastic problem (1) as an abstract form, with the nonlinear operator ∶ R n → R n , the stochastic solution u (t, ) ∈ R n×n t , the stochastic force F (t, ) ∈ R n×n t and n is the number of degrees of freedom of the spatial discretization, n t is the number of total time steps. In practical numerical implementation, the matrix/vector implementations of (u (t, ) , ) and F (t, ) are calculated by using (6) and (7), respectively. To solve the stochastic solution of (20), we extend the Newton (or called Newton-Raphson) method for solving deterministic nonlinear equations to stochastic cases. A stochastic Newton method is used to linearize (20) as for i = 1, … , n t , where the stochastic tangent stiffness matrix K T (u k−1 (t i , ) , ) ∈ R n×n and the stochastic force vector F k (u k−1 (t i , ) , ) ∈ R n are given by and the matrix assembly of K T (u k−1 (t i , ) , ) (or (u) ∕ u) is implemented by using the tangent tensor  T in (19). In this way, the stochastic solution u k (t i , ) of the kth iteration at the time step t i in (21) is approximated as where u k−1 (t i , ) ∈ R n is the stochastic solution of the previous iteration that has been determined and Δu k (t i , ) ∈ R n is the kth stochastic increment that needs to be solved.

Approximation of stochastic solutions
In practice, (24) is only of theoretical significance since Δu k (t i , ) is not as easy to be solved as deterministic cases. On the basis of (24), we develop a more detailed approximation for the stochastic solution u (t i , ) of each time step t i , which is given by whereũ (t i−1 , ) is the stochastic solution of the time step t i−1 and has been determined, is the corresponding matrix, k is the number of approximate terms. Equation (25) is a kind of intrusive approximation, but we approximate the random variables j (t i , ) using random samples in practice, which is a non-intrusive way. Combining these two points, (25) is considered as a weakly intrusive approach and we will explain this point in detail in the subsequent section. In fact, the solution u k (t i , ) can be considered as a random field. Equation (25) is very similar to the classical series expansions of the random field u k (t i , ), for example, the Karhunen-Loève expansion and the Polynomial Chaos expansion 11 that decompose the random fields into stochastic and deterministic spaces. Different from the known probability properties such as mean values and covariance functions of the classical random fields, the constraints of u k (t i , ) are only from the stochastic governing Equation (20). In this way, all quantities } k j=1 cannot be known a priori. The random variables } k j=1 may no longer be uncorrelated as they are in the Karhunen-Loève expansion. We need to address this issue and develop dedicated numerical methods for the purpose. The convergence and accuracy of (25) may benefit from these classical expansions, but are not well understood currently. Further investigation is still required. Several explorations for solving linear static stochastic problems, linear dynamic stochastic problems and the problems defined on random domains using (25)-like expansions can be found in our previous work. 30,31,33 We also remark that the second part-like expansion in (25) ) has been widely used in the context of proper generalized decomposition methods. [34][35][36] However, these classical expansions are not applicable to the stochastic elastoplastic problem. In this article, we develop an efficient and accurate numerical algorithm to solve stochastic elastoplastic problems by using the weakly intrusive approximation (25) and a dedicated iteration.
Combining (24) and (25) we rewrite (25) as an incremental format similar to (24). To this end, the previous approximation u k−1 (t i , ) and the stochastic increment Δu k (t i , ) in (24) are given by where the first k − 1 couples are assumed to be known. In this way, we can still solve the unknown stochastic increment by taking advantage of the stochastic Newton method mentioned in the previous section and the stochastic increment is given by (27) and (28) into (21) we have An iterative algorithm for solving the couple { k (t i , ) , d k (t i )} will be discussed in the next section.

Iterative algorithm for stochastic solutions
It is not easy to calculate k (t i , ) and d k (t i ) in (29) at once. We adopt an alternating iterative method to solve k (t i , ) and d k (t i ) one by one. Specifically, we are to fix one of them and solve the other, and then the fixed one is updated according to the solution of the other. In detail, if the random variable k (t i , ) has been known (or given an initial value), the stochastic Galerkin method 11,12 is used to transform the original stochastic problem (29) into a deterministic finite element equation where the deterministic matrix K * T,k (t i ) ∈ R n×n and the deterministic vector F * k (t i ) ∈ R n are calculated by where E {⋅} represents the expectation operator. It is noted that the above equations are the approximations of the proba- where  ( ) is the probability measurement of the random input. In this way, the probability integrals can be achieved cheaply, but theoretically, high-accuracy approximation requires a large number of random samples. In numerical examples, we will show that the sample size requirement can be greatly relaxed, that is, only a small number of random samples are used to perform (31) and (32), but high-accuracy stochastic solutions are still obtained. Equation (30) is a linear finite element equation and existing linear solvers 4 can be adopted to efficiently solve the vector d k (t i ). To speed up the convergence in practical computations, we orthogonalize the vector d k (t i ) by Gram-Schmidt orthogonalization, which corresponds to where are orthogonal vectors that have been determined via the first k − 1 iterations and they meet According to the deterministic vector d k (t i ) obtained by (30), the random variable k (t i , ) is updated by a classical Galerkin procedure where the random variables a k (t i , ) ∈ R and b k (t i , ) ∈ R are given by A common method for solving (34) is known as the PC-based methods. 11,37 In this method, the random variable k (t i , ) is expanded by PC bases and the expanded coefficients are calculated via an augmented deterministic equation. It alleviates the computational burden as the augmented deterministic equation is not coupled with degrees of freedom discretized by the finite element method, but it does not overcome the curse of dimensionality arising in high-dimensional stochastic problems. To avoid this issue, a non-intrusive method is adopted to calculate random samples of the random variable

Improved scheme of stochastic solutions
We recall the approximation (28) and the stochastic increment is approximated by a stochastic-deterministic decoupling decomposition Δu , which has low accuracy in many cases. For the explanation of this point, we consider the sample descriptions Δu (28) is then written as which is considered as a rank-1 approximation of the matrix Δu k ( t i ,̂), thus it is less accurate for many problems. To avoid this problem, we let the known matrix D (t i ) in (26) be a set of reduced bases and recalculate the random vector (t i , ). To solve the random vector (t i , ), we still adopt the stochastic Newton method for the linearization and the Galerkin method for building a reduced-order linearized stochastic finite element equation, where the vector (j) (t i , ) ∈ R k is the stochastic solution of the jth iteration and the vector Δ (j) (t i , ) ∈ R k is the jth stochastic increment. To solve the increment Δ (j) (t i , ), we use a sampling method to solve (39) for all sample realizations, which is achieved by solving the following equatioñ for q = 1, … , n s , where the reduced-order matrixk T ( (j) ( t i , (q) )) ∈ R k×k and the reduced-order vector Further, the random vector is updated via The computational effort of (40) is very low since its size k is usually not too large. For all sample realizations { Δ (j) ( t i , (q) )} n s q=1 , total computational costs of n s solutions of (40) are still low. It is noted that compared to intrusive methods, the non-intrusive approximation of the random vector (t i , ) more easily and accurately captures the nodes in elastic/plastic states. Also, another advantage of the proposed method is that the initial random samples of random sources are propagated from beginning to the end, that is, the sample realization (q) is propagated to the final stochas- , which provides an easy-to-use and highly accurate way for uncertainty quantification.

Computational aspects
In this section, we propose efficient numerical formats to accelerate the numerical computation of the proposed method. Computational costs of the proposed method are concentrated in the assembly of the stochastic tangent stiffness matrix (29), the expectation calculations in (31) and (32), the vector-stochastic-matrix multiplications in (35) and (36), and the matrix-stochastic-matrix multiplications in (41) and (42). To speed up these computations, we follow the idea for numerical implementation in Reference 38 and write the stochastic tangent stiffness matrix where B ∈ R n I ×n is a deterministic sparse matrix related to the strain-displacement operator, n I is the total number of integration points, M T ( ) ∈ R n I ×n I is a block diagonal sparse matrix related to the stochastic elastoplastic constitutive model. We decompose the stochastic matrix M T ( ) into two parts where the stochastic matrix M e ( ) ∈ R n I ×n I is related to the elastic constitutive part (i.e., linear problem) of the material and the stochastic matrix M T,p (u k−1 (t i , ) , ) ∈ R n I ×n I is related to the elastoplastic constitutive part (i.e., nonlinear problem). The stochastic matrix M T,p (u k−1 (t i , ) , ) is much sparser than the stochastic matrix M T (u k−1 (t i , ) , ) since the nodes in plastic states are only a part of the whole domain, or even a small part in most cases. Thus, the stochastic matrix K T (u k−1 (t i , ) , ) can be rewritten as where the stochastic elastic stiffness matrix K e ( ) ∈ R n×n corresponding to M e ( ) and the stochastic plastic tangent stiffness matrix K T,p (u k−1 (t i , ) , ) ∈ R n×n corresponding to M T,p (u k−1 (t i , ) , ) are Similarly, the stochastic force vector can be written as It is noted that the matrix K e ( ) can be precomputed via M e ( ) and does not need to be updated during the iterations, while the matrix 1 (t i , ) , ) (i.e., the stochastic matrix corresponding to plastic points), where n Ip ( ) is a random variable representing the number of integration points of elements in plastic states. In the numerical implementation, we compute and store M nonzero ̂) ,̂)), which only requires little computational effort and storage memory. In this way, we can efficiently implement the proposed method. Specifically, (31) and (32) are calculated via whereÊ {⋅} represents the operator for estimating expectations using random samples, for instance, we havê Further, the random sample vectors of the random variables a k (t i , ) and b k (t i , ) in (35) and (36) are calculated by where the deterministic vector Bd k (t i ) ∈ R n I is only calculated at once. Also, introducing a deterministic matrix B D (t i ) = BD (t i ) ∈ R n I ×k we calculate (41) and (42) bỹ

Algorithm scheme
The proposed nonlinear stochastic finite element method for solving stochastic elastoplastic problems is summarized in Algorithm 1, which consists of three-loop procedures. The outer loop is from step 3 to step 26, which is used to compute the stochastic solution u (t i , ) at each time step t i and the looping number of which is the same as total time steps n t . The middle loop from step 5 to step 20 is used to solve the couples (27), and the inner loop from step 8 to step 16 is used to solve each couple { k (t i , ) , d k (t i )}. Before the outer loop, we precompute the stochastic matrices M e (̂) ∈ R n I ×n I ×n s,1 , M e (̂ * ) ∈ R n I ×n I ×n s,2 and the deterministic matrix B ∈ R n I ×n I in step 1, where the sample size n s,1 is used to solve the deterministic vector { d j (t i ) } k j=1 in the middle and inner loops and the sample size n s,2 is adopted to solve the final stochastic solution in the outer loop. Significant computational effort is saved by setting n s,1 ≪ n s,2 , but good accuracy is still achieved. The performance of this strategy will be illustrated in numerical examples. In step 2, stochastic solutionsũ (1) ( t 0 ,̂) andũ while u,k ≥ u do 6: Calculate the sparse stochastic matrix M T,p Initialize the random samples (0) ∈ R n s,1 8: while d,j ≥ d do 9: Assemble the matrix K * T,k (t i ) and the vector F * k (t i ) by (49) and (50) 10: Solve the deterministic vector d Orthogonalize d (j) Calculate the random sample vectors a k (52) and (53) 13: Solve the random sample vector Compute the stochastic solutionũ (2) ( t i ← t i+1 26: end for that any nonzero vectors of size n s,1 can be used as the initialization and it almost has no influence on the computational accuracy and efficiency of the proposed method. With the known random sample vector, the deterministic vector d are recalculated by solving n s,1 and n s,2 k-dimensional nonlinear stochastic finite element equations using Algorithm 2, which simply extends the classical Newton iteration to the reduced-order nonlinear stochastic finite element equations but has better convergence and lower computational effort. Corresponding stochastic solutionsũ (1) ( are thus updated. The stochastic solutionũ (1) ( t i ,̂) is only used to the next middle and inner loops andũ is considered as the final stochastic solution of the time step t i .

Algorithm 2.
Algorithm of the reduced-order Newton iteration 1: for q = 1, … , n s do 2: Initialize the solution ( (q) ) = 0 ∈ R k 3: while NR,k ≥ NR do 4: Assemble the reduced-order tangent matrixk T ( (j) (t i , ) ) ∈ R k×k and the vectorf (k) ( (j) (t i , ) ) ∈ R k by (54) and (55) 5: Solve the increment Δ ( (q) ) ∈ R k by (40) 6: Update the solution (k+1) ( (q) ) = (k) ( (q) ) + Δ (k) ( (q) ) 7: Compute the iterative error NR,k 8: k ← k + 1 9: end while 10: q ← q + 1 11: end for Two iterative criteria in Algorithm 1 are used to check the convergence, that is, d,j in step 14 and u,k in step 18. The locally iterative error d,j is defined as which measures the difference between the vectors d Similarly, the globally iterative error u,k is defined as that are calculated as following. 31 To this end, we calculate the autocorrelation function of the random vector (t i , ) by which is decomposed into by the eigen decomposition, where Q (t i ) ∈ R k×k is an orthonormal matrix and Z (t i ) is a diagonal matrix with elements on the diagonal in descending order. A new random vector is given bỹ( ] T ∈ R k and its autocorrelation function is C̃̃(t i ) = Z (t i ). Substituting new random variables into (57) the iterative error u,k becomes as where Tr (⋅) is the trace operator and Z k (t i ) is the kth diagonal element of the matrix Z (t i ). In this way, the iterative error u,k keeps decreasing as the retained item k increases. More details of this improved error estimation can also be found in Reference 31. Also, the iterative error NR,k in step 3 in Algorithm 2 is defined as which measures the contribution of the kth increment Δ (k) ( (q) ) to the solution (k) ( (q) ) and the iteration stops when the increment contributes little to the final solution.
Let us revisit the weak intrusiveness of the proposed method. On one hand, the matrix K * T,k (t i ) and the vector F * k (t i ) assembled in step 9 in Algorithm 1 only involve the sample matrices M e (̂) , ̂) ,̂) and the sample vec- (49) and (50), which can be assembled using existing codes in non-intrusive ways (a slight modification taking advantage of sparsity to save memory and accelerate computation). In step 10, existing numerical solvers are used to solve (30), where the matrix K * T,k (t i ) keeps the same size and matrix properties as the original stochastic matrix K T (u k−1 (t i , ) , ). On the other hand, the recalculated Equation (40) only requires few effort for assembling the reduced-order matrix and vector in (41) and (42) and its solution by Algorithm 2 is fully non-intrusive. In these senses, we implement the intrusive approximation (25) of the stochastic solution only in weakly intrusive ways and most existing codes can be reused as-is or slightly modified.

NUMERICAL EXAMPLES
In this section, we test the proposed nonlinear stochastic finite element algorithm by a two-dimensional stochastic elasto-

Problem setting
We consider a plate with a hole shown in Figure 1A, where the length l = 2 m and the width w = 1 m. The radius of the circle is r = 0.2 m. Vertical forces f (x, y) are applied on two sides. We only consider a quarter plate shown in Figure 1B for the analysis due to the symmetry. It is noted that we need to assume that the random material properties have the same symmetry as the geometry. Otherwise, we can only perform the analysis on the whole domain instead of the quarter. By means of finite element discretization, n p = 292 nodes and n e = 514 linear triangular elements are generated, and n = 584 degrees of freedom are involved in total. The stochastic displacement in the x direction u x (t, ) = 0 holds on the boundary Γ D1 and the stochastic displacement in the y direction u y (t, ) = 0 holds on the boundary Γ D2 . We only consider a deterministic hardening parameter ( ) = 1 × 10 4 Pa in (11), but it is not a limitation of the proposed method. A random variable/field ( ) is easy to be embedded into the proposed computational framework. The stochastic bulk modulus ( ) and the stochastic shear modulus ( ) in (10) are given by where the Poisson rate = 0.29, the Young's modulus E (x, y, ) is a Gaussian random field with the mean value E 0 (x, y) = 2.07 × 10 5 Pa and the covariance function where the standard deviation Cov = 0.1E 0 (x, y) and the correlation lengths l x = l y = 1 m. By using the Karhunen-Loève (KL) expansion, 39,40 the random field E (x, y, ) is approximated as where r E is the number of the truncation, { j ( ) } r E j=1 are mutually independent standard Gaussian random variables, and { E,j , E j (x, y) } r E j=1 are eigenvalues and eigenfunctions of the covariance function Cov (x 1 , y 1 ; x 2 , y 2 ). They are solved by the following homogeneous Fredholm integral equation of the second kind, Similarly, the yield stress Y (x, y, ) is a Gaussian random field with the mean value Y,0 (x, y) = 250 Pa and the covariance function given in (63). It is also approximated as where the standard deviation Cov = 0.1 Y,0 (x, y) and are mutually independent standard Gaussian random variables and they are independent of the random variables . We let the number of the truncation r Y = r E = 10 and the stochastic dimensionalities are 20 in total. Studies on the truncation error can be found in Reference 25, which is beyond the scope of this article. According to (65), Y ,i = E,i and Y,i (x, y) = E i (x, y), i = 1, … , r E hold. The samples (i) such that min x,y∈Ω E ( x, y, (i) ) ≤ 1 × 10 −3 and min x,y∈Ω Y ( x, y, (i) ) ≤ 1 × 10 −3 are dropped out in practical implementation. Furthermore, in many applications, the modified exponential covariance kernel is suggested to simulate random fields due to its good properties. 41,42 We only consider the classical exponential covariance kernel in this article. It is easy to consider the modified exponential covariance kernel in the proposed SFEM by a slight modification of (63).
According to the KL expansion (64), the stochastic matrix M e ( ) in (45) has a separated form The elastic stochastic stiffness matrix K e ( ) is thus calculated by where the deterministic matrices K e,i = B T M e,i B ∈ R n×n , i = 0, 1, … , r E . The matrices in (49), (52), and (54) are then computed via In this way, we can deal with the stochastic elastic part of the nonlinear problem by precomputing and storing the matrices B and { K e,i } r E i=0 . Only the stochastic matrix M p ( ) needs to be assembled in each iteration, which saves a lot of computational effort and storage memory.

Example 1: Subjected to a monotonic load
In this case, we consider the force f (x, y) = 200N, which is loaded monotonically in n t = 11 time steps {t 0 , t 1 , … , t 10 }. For each time step, the iterative errors u,k are shown in Figure 2, which demonstrates that the proposed method has good convergence. For all time steps, k = 6 to 12 retained terms

F I G U R E 3 The first six vectors
in the y direction at the time step t 1 (the first line) and the time step t 10 (the second line).  Figure 1B) obtained by the proposed nonlinear SFEM and MCS are compared in Figure 4A, which indicates that the PDFs obtained by the proposed method are very consistent with that of MCS. To quantify the approximate errors of PDFs, the absolute errors between the proposed SFEM and MCS are given in Figure 4B. Since the values of tails of the PDFs are very close to zero, relative errors are not suitable to describe the approximate errors for whole PDFs. The relative error around the maximum PDF (the corresponding displacement is about 7.0 × 10 −4 ) of the point A can be roughly estimated to be 1%, and that of the point B (the corresponding displacement is about 3.2 × 10 −4 ) is less than 1%, which indicates the high accuracy of the proposed method. We highlight that the proposed method can also be used as a stochastic reduced-order method. To show this point, ten sample realizations of displacement-load curves of the points A and B obtained by the reduced basis { d j (t i ) } j=1 and MC simulations (i.e., full-order equations) are depicted in Figure 5A,C, respectively, and the corresponding errors relative to MCS are shown in Figure 5B,D. The maximum relative errors of both point A and point B are less than 0.01%, which demonstrates the high-accuracy approximation of the reduced-order equations.  1 (t i , ) , ) and the finite element Equation (30) of the fine mesh are larger than that of the coarse mesh, higher computational effort is required for the fine mesh.

Example 2: Subjected to a cyclic load
In this case we consider a cyclic load shown in Figure 7 and n t = 41 time steps {t 0 , t 1 , … , t 40 } are adopted for the time discretization. The numbers of retained terms at different time steps {t i } 40 i=1 are listed in Figure 8. It is seen that more retained  terms are required for the time steps close to t 10 and t 30 . Since the time steps t 10 and t 30 correspond to the maximum and the minimum loads, more terms are needed to approximate the strongly nonlinear stochastic solutions. As comparisons, we need less retained terms to approximate the stochastic solution of the step times close to t 20 Figure 9 are only used to capture the stochastic displacements caused by the stochastic residual plastic strains, and its number is only slightly more than the number of vectors for approximating stochastic solutions of elastic states.
The PDFs of the displacements of the points A and B in the y direction at the time step t 40 are shown in Figure 10A and the corresponding absolute errors are shown in Figure 10B. The relative errors around the magnitudes of PDFs (the corresponding displacements are about −4.5 × 10 −5 for the point A and −6.3 × 10 −5 for the point B) are less than 2.5% and 2.0%, respectively, which indicates that the proposed SFEM still has high accuracy for the case under cyclic loads. But the errors are slightly increased than that of the case under monotonic loads illustrated in Example 5.2. Furthermore, ten sample realizations of displacement-load curves of the points A and B obtained by the proposed nonlinear SFEM and MC simulations are shown in Figure 11A,C, respectively. The corresponding errors relative to MCS are shown in Figure 11B,D, which indicates that relative errors of stochastic displacements of the points A and B are very small and the proposed SFEM can still be used as a stochastic reduced-order approach under complex loading. It is also seen that most relative errors are still less than 0.01%. The relative errors are a bit larger for some steps corresponding to the time steps 20 to 30 due to the stochastic displacements close to zero. The computational times of the proposed nonlinear SFEM and MC simulations are 427.88 s (including the solving time 88.16 s and the updating time 339.72 s) and 6414.12 s, respectively, which demonstrates that the proposed method is still cheap enough. Also, compared to the computational cost for this example and that for Example 5.2, it is found that the computational cost almost proportionally increases with the number of total time steps. Thus, reducing the number of total time steps with an adaptive time-stepping method will greatly reduce the computational cost, which is beyond the scope of this article and will be investigated in future research.

Example 3: High-dimensional stochastic case
In this example, we consider a high-dimensional stochastic case by letting the number of the truncation r E = r Y = 50 in (64) and the total stochastic dimensionalities are 100 in this case. The sample size in step 2 in Algorithm 1 are reset  as n s,1 = 1000. The monotonic load in Example 5.2 and 11 time steps are adopted. The iterative errors u,k of each time step {t i } 10 i=1 are shown in Figure 12 and k = 8 to 14 retained terms } k j=1 are required for the specified precision. Compared to the iterative errors shown in Figure 2, the convergence rate is lower in this case due to the large uncertainties caused by high-dimensional random variables. The comparison between the number of retained terms of each time step for the stochastic dimensionalities 20 (i.e., the case in Example 5.2) and 100 is depicted in Figure 13. The number of retained terms of each time step in this case is greater than that of the stochastic dimensionalities 20 and 109 deterministic vectors d j (t i ) ∈ R 584 and 109 random sample vectors j ( t i ,̂) ∈ R 10 4 are involved in total. The solving and updating times of the proposed method and the MCS time are 49.92, 137.62, and 1726.35 s, respectively, which is a little bit time-consuming since more retained terms } are required to be resolved. The PDFs of the stochastic displacements of the points A and B in the y direction at the time step t 10 are shown in Figure 14A and the corresponding absolute errors are seen from Figure 14B. Comparing Figures 4B and 14B, larger absolute errors are obtained in high-dimensional stochastic cases. The PDFs shown in Figure 14A are very close to that of the low-dimensional case considered in Figure 4A since the last few dozens of random terms expanded by KL expansion (64) maybe have less contributions to the quantities of interest. It is noted that the proposed SFEM solves both low-and high-dimensional stochastic problems in a unified computational frame and no additional modifications for    high-dimensional stochastic problems are performed. Thus, we only need to execute the stochastic analysis under 100 stochastic dimensionalities directly in this case. If there are larger uncertainties caused by high-dimensional expansions, only more additional retained terms are required to capture the large variability of the stochastic solutions.

Example 4: Random sources with large variability
In this example, we consider that the stochastic Young's modulus E (x, y, ) and the stochastic yield stress Y (x, y, ) have larger variability by letting the standard deviations as 0.2E 0 (x, y) and 0.2 Y,0 (x, y), that is, 0.2 times the mean values, hereinafter which is simply notated as Cov = 0.2□. The number of the truncation r E = r Y = 10 in (64) and the total stochastic dimensionalities are 20. The monotonic load in Example 5.2 and 11 time steps {t 0 , t 1 , … , t 10 } are still adopted in this case. The comparison between the number of retained terms of each time step for the cases Cov = 0.1□ (i.e., the case in Example 5.2) and Cov = 0.2□ is depicted in Figure 15. It is seen that the number of retained terms of the case Cov = 0.2□ are greater than that of the case Cov = 0.1□ for all time steps, and 139 terms  Since more terms } are solved in this case, the computational cost is little higher than that of the case Cov = 0.1□. The PDFs of the displacements of the points A and B in the y direction at the time step t 10 are shown in Figure 16A and the corresponding absolute errors are shown in Figure 16B. Magnitudes of absolute errors are very small compared to the magnitudes of their PDFs and the proposed SFEM still has good accuracy in this case. Also, the proposed method can well capture the long-tailed characteristics of the PDFs, which is of great help for uncertainty quantification in many stochastic problems. For instance, a high-accuracy approximation of PDF tails is critical in structural reliability analysis. Applying the proposed SFEM to uncertainty quantification for stochastic elastoplastic problems is beyond the scope of this article and will be carried out in a follow-up study.

CONCLUSIONS AND DISCUSSIONS
This article develops an accurate and efficient nonlinear stochastic finite element method to solve stochastic elastoplastic problems. A stochastic Newton method is adopted to linearize the nonlinear stochastic finite element equation at each time step. A weakly intrusive approximation of the stochastic solution and a dedicated iteration are used to solve the corresponding linearized stochastic finite element equations. In this way, the original stochastic problem is transformed into a series of deterministic linear finite element equations and the corresponding one-dimensional stochastic algebraic equations. Existing FEM solvers can be used to solve the deterministic linear finite element equations without any modification. All stochastic analyses are performed by the one-dimensional stochastic algebraic equations that are efficiently solved by a proposed non-intrusive method. The curse of dimensionality of high-dimensional stochastic elastoplastic problems is thus circumvented with great success, which has been illustrated via a numerical case of up to 100 stochastic dimensionalities. However, there are still several intriguing issues that require to be investigated further. Only stochastic Newton iteration is used to deal with the stochastic nonlinearity in this article. Other nonlinear iterations, for example, stochastic quasi-Newton method, can also be developed for this purpose. Also, it is attractive to introduce adaptive time-stepping strategies to reduce the number of total time steps, which will save a lot of computational effort of the proposed method.

ACKNOWLEDGMENTS
The authors are grateful to the Alexander von Humboldt Foundation and the International Research Training Group 2657 (IRTG 2657) funded by the German Research Foundation (DFG) (Grant number 433082294). Open Access funding enabled and organized by Projekt DEAL.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.