Effect of Element Size in Random Finite Element Analysis for Effective Young ’ s Modulus

In random finite element analysis (RFEA), continuous random fields must be discretized. The critical element size to achieve acceptable accuracy in effective Young’s modulus for an elementary soil mass is investigated. It is observed that the discrepancy between the continuous and discretized solutions is governed by the discretization strategy (element-level averaging versus midpoint), spatial variability pattern, and the adopted autocorrelation function. With the element-level averaging strategy, RFEA with element size less than (scale of fluctuation)/5 will not induce significant discrepancy from the continuous solution. Moreover, the element-level averaging strategy is more effective than the midpoint strategy.


Introduction
Spatial variability of soil parameters is commonly encountered in naturally occurring materials such as soils and rocks.Behaviors of foundations and slopes in spatially variable soils can be simulated by random finite element analysis (RFEA) or random finite difference analysis [1][2][3][4][5][6][7][8][9][10][11].An inevitable approximation in RFEA is that the continuous random field must be discretized, as the soil parameter within each element must be homogeneous in common finite element packages.The RFEA solution with such an approximation is in general different from the continuous solution.
The continuous random field can be discretized in several possible ways.A well-known choice in the literature is to adopt the local spatial average of the random field over each element [12,13].Fenton and Vanmarcke [14] developed an algorithm called the local averaging subdivision (LAS) to implement this local spatial averaging procedure, whereas Jha and Ching [15] showed that a Fourier series method can achieve the same purpose.Alternatively, the midpoint strategy [16] is to take the random field value at the centroid of each element.There are also strategies based on shape functions [17], series expansions [18], and optimal linear estimation [19].Li and Kiureghian [19] conducted a review of these methods and compared their accuracy in representing the continuous random field, quantified by the variance for the maximum error between the discretized and continuous fields.More recently, Ching and Phoon [20] and Huang and Griffiths [21] considered the shear strength of a twodimensional soil column and investigated the discrepancy between the discretized and continuous solutions and its relationship with RFEA mesh size.Their conclusions are quantitatively different.In order to keep the discretization error acceptable, [20] recommended a mesh size that is an order of magnitude smaller than that recommended in [21].The reason for this significant difference is under investigation, but it should have something to do with the definition of "acceptable accuracy." In contrast to [20,21], the current study focuses on the discrepancy between the discretized and continuous solutions for "effective Young's modulus."Young's modulus of a soil mass is more relevant than its shear strength for designs governed by the serviceability limit state.In fact, the need to pay some attention to the serviceability limit state is reflected by the renaming of ISSMGE TC23 "Limit State Design in Geotechnical Engineering" to TC205 "Safety and Serviceability in Geotechnical Design" around 5 years earlier.
In this study, effective Young's modulus ( eff ) is defined as overall Young's modulus of a spatially variable elementary soil mass subjected to one-dimensional compression (similar to that in an oedometer cell).Loading is applied through a prescribed uniform displacement at one boundary for the elementary soil mass.
The discrepancy between the discretized and continuous solutions for  eff is inevitable because  eff of an element is, for example, neither the spatial average over the element nor the random field value at the centroid.The normalized discrepancy between the discretized and continuous solutions for  eff is quantified as follows: [(statistics for discretized  eff ) -(statistics for continuous  eff )]/(statistics for continuous  eff ).Only two statistics, namely, mean and standard deviation, are considered.Practical guidelines on minimum element size needed to achieve an acceptable normalized discrepancy less than 0.01 and 0.05 are presented based on extensive RFEA.The discrepancy between the discretized and continuous solutions for  eff can be minimized by adopting very small elements in RFEA, because the difference between  eff and the spatial average (or the centroid value) diminishes as the element size decreases.Nonetheless, it may not be realistic to adopt very small elements due to limited computation resources.
The focus of this paper is therefore to clarify the following: (1) what the relationship between the discrepancy and element size is, (2) whether this relationship would change with factors such as the degree of random field anisotropy and loading direction, and (3) how the conclusion compares with that for shear strength, for example, [20,21].It is recognized that the mechanical properties of a rock mass change depending on the size of the specimen due to the presence of discontinuities.The random fields studied in this paper do not produce spatial variations as localized and contrasting to behave in the manner of discontinuities.Numerical studies with varying specimen sizes further confirm that scale effect is not present in the results presented below.

Random Finite Element Model
2.1.Random Field Model.Spatial variabilities of soil properties are usually modeled by random fields [12].Among random field models, stationary random fields are widely used due to their simplicity and possibly the only practical version that can be characterized from limited data [22].Consider a spatially variable soil mass in a two-dimensional (2D) plane strain condition.The size of the domain is  ×  = 12.8 × 12.8 (Figure 1).No unit is specified for the dimension L because all dimensions will be scaled by the scale of fluctuation (SOF).There is a potential for scale effect for random fields producing features similar to discontinuities (extremely localized and contrasting features), but this paper only examines isotropic and layered anisotropic random fields.Young's modulus within the soil mass, denoted by (, ), is simulated as a stationary lognormal random field with inherent mean =  = 100 and inherent coefficient of variation (COV) =  = 0.5, where (, ) are the horizontal and vertical coordinates.Therefore, ln[(, )] is a stationary normal random field.No unit will be specified for  because all statistical properties of the soil mass will be later scaled by .To define the correlation structure in ln[(, )] between two locations with horizontal distance = Δ and vertical distance = Δ, the single exponential (SExp) autocorrelation model is considered [12,23]: where   and   are the horizontal and vertical SOFs, respectively, for the ln[(, )] random field.Another popular autocorrelation model is the squared exponential (QExp) model: It is clear that the correlation decreases as  and  increase.Because of the manner in which most natural soils are deposited, soil properties are expected to be strongly correlated within a small interval and weakly correlated when measurement points are separated sufficiently far apart.
The first author has developed the Fourier series method (FSM) for simulating stationary normal random fields [15,24].A 2D stationary lognormal random field (, ) can be simulated by taking the exponential of a 2D stationary normal random field ln[(, )] with mean where Re[⋅] denotes the real part of the enclosed complex number;   and   are independent zero-mean normal random variables with variance  2  given by [15] where   = /  and   = /  .Besides simulating the point process of the normal random field ln[(, )], the FSM is also able to directly simulate the arithmetic averages of the normal random field ln[(, )] over rectangular regions (cells) in 2D [15].In this study, the cells are the small elements in the RFEA.The "local average" for each small element is taken to be the exponential of the arithmetic average for ln[(, )] over that element.By doing this, the local average for (, ) is actually a geometric average.This local averaging will be later referred to as "element-level averaging" (ELA).Poisson's ratio (]) is set to be a constant (] = 0.3) because the impact of the spatial variability of Poisson's ratio is insignificant [2,25,26].

Finite Element
Model.The 12.8 × 12.8 square domain (an elementary soil mass) is modeled by the RFEA mesh shown in Figure 1.The commercial code ABAQUS is adopted.Fournode plane strain elements with full integration (CPE4) are used.Each element follows isotropic elasticity with E = ELA value and with ] = 0.3.For each realization of the E random field, two RFEAs that simulate an overall 1D compression in the x or z directions are conducted.Let us consider the x direction as an example.As shown in Figure 1, the nodes on Side 1 are constrained in a "tied freedom" manner [27] so that all nodes on Side 1 are subjected to a uniform xdisplacement of 0.1 to the left.Sides 2 to 4 are constrained by rollers (e.g., for Face 2, z-displacement is not allowed).The predominant strain in the soil mass should be in the x direction.The other strain components,   and   , may not be zero in the soil mass due to spatial variability, but they are likely to be minor.The boundary condition described above is similar to that imposed in an oedometer test.Moreover, the tied freedom approach allows an exact back-calculation of effective elastic parameters [27].The effective Young modulus in the x direction, denoted by  ,eff , can be deduced from the response of the elementary soil mass.The stress at Side 1 is not uniform because a uniform displacement boundary is prescribed.The "overall"   is equal to the arithmetic average   over Side 1.Because the overall strains are   = 0.1/12.8and   =   = 0 (in the global sense, not strain at a point), it can be shown that effective (overall) Young's modulus for the x direction is It should be emphasized that  ,eff is related to the overall 1D compression response of the soil mass in the x direction under ] = 0.3.A similar displacement-controlled loading is applied in the z direction over the same realization of the E random field.Therefore, each E random field realization produces a set of ( ,eff ,  ,eff ) samples.

Random Finite Element Analyses with Variable Element Sizes
Discrepancy between the continuous solution of  eff and any discretized solution will be present.Let the continuous solution be denoted by  eff and let the RFEA solution be denoted by  eff RFEA .In this section, detailed analyses will be carried out to understand the magnitude of the normalized discrepancy between  eff and  eff RFEA and its relationship with the element size.The RFEA with variable element sizes are taken: eight element sizes corresponding to 0.1, 0.2, 0.4, 0.8, 1.6, . .., and 12.8, respectively, denoted by Size 1, Size 2, . .., and Size 8. Size 1 mesh is the finest with element size = 0.1 × 0.1.There are therefore 128 × 128 = 16,384 elements.Size 8 mesh is the coarsest with element size = 12.8 × 12.8.There is only a single element.The RFEA starts from Size 1 mesh.Size 1 mesh is then subjected to displacement-controlled compression loadings in x and z directions to obtain  ,eff RFEA and  ,eff RFEA .Size 2 mesh uses Size 1 mesh but conducts geometric averaging over 2 × 2 = 4 elements in Size 1 to obtain the E values for Size 2 (see Figure 2).Size 2 mesh is then subjected to displacementcontrolled compression loadings to obtain  ,eff RFEA and  ,eff RFEA .The same procedure is repeated to progressively obtain the meshes for Size The normalized discrepancy is defined as [(statistics of   ) − (statistics of  1 )]/(statistics of  1 ): it quantifies the discrepancy between the statistics of   and the statistics of  1 .Note that  1 serves as the comparison baseline because it is assumed to be close to the continuous solution  eff .Figure 3 shows the logic block diagram for the main task.The recommendations about the mesh sizes will be given based on whether the absolute value of the normalized discrepancy is less than a prescribed threshold (0.01 or 0.05).Normalized discrepancy:

Isotropic Cases.
Normalized discrepancy:   4(a).The rightmost "◻" corresponds to Size 1 mesh [/(element size) = 1000], whereas the leftmost one corresponds to Size 8 mesh [/(element size) = 0.78].For isotropic cases, the results for    and    are the same.Therefore, only the results for    are presented in Figure 4.
For isotropic cases, the discrepancy in the sample mean is always positive (i.e.,    sample mean greater than   1 sample mean).This implies that the use of a coarse mesh in place of a fine mesh will overpredict  eff (unconservative).
The discrepancy in the sample STD is also positive (i.e.,    sample STD greater than   1 sample STD).However, both discrepancies are quite small: their absolute values are less than 0.05 even when /(element size) is very small (or element size is very large).The small absolute normalized discrepancies are expected even for coarse meshes, because [2,25,26] confirmed that, for isotropic cases,  eff of a domain can be well approximated by its geometric average.Recall that in the RFEA the element-level average (ELA) is geometric average.This means that even for the coarsest mesh (element size = 12.8), the discretization is a satisfactory approximation.Therefore, although the use of a coarse mesh in isotropic cases will overpredict  eff , the degree of overprediction is quite slight.According to Figure 4, if the acceptable absolute normalized discrepancy is 0.01, /(element size) needs to be greater than 5.This means that the element size needs to be smaller than /5 in order to have |normalized discrepancy| < 0.01.This ratio of 5 is called the critical ratio (  ) for an acceptable error = 0.01.

Layered Cases. For layered cases, Figures 5(a) and 5(b)
show the normalized discrepancies in sample mean and STD for    , effective Young's modulus when the elementary soil mass is subjected to x direction displacement-controlled loading (parallel to layers).In contrast to isotropic cases, the discrepancy in the sample mean is always negative (i.e.,    sample mean smaller than   1 sample mean), and the absolute normalized discrepancies are larger than those for isotropic cases (Figure 4(a)).The larger discrepancies are expected, because [26] confirmed that for layered cases,  ,eff can be well approximated as its arithmetic average.However, the ELA adopted in this study is geometric average, not arithmetic average.If the acceptable error is |normalized discrepancy| < 0.01, /(element size) needs to be greater than 5 (  = 5).
Figures 5(c) and 5(d) show the normalized discrepancies in sample mean and STD for    (displacement direction perpendicular to layers).The absolute normalized discrepancies are significantly larger than those for isotropic cases (Figure 4) and also larger than those for layered  ,eff (Figures 5(a) and 5(b)).The discrepancy in the sample mean is always positive (i.e.,    sample mean larger than   1 sample mean).Ching et al. [26] confirmed that for layered cases  ,eff can be well approximated by its harmonic average.However, the ELA adopted in this study is geometric average, not harmonic average.If the acceptable error is |normalized discrepancy| < 0.01, /(element size) needs to be greater than 7.8 (  = 7.8).
Column (a) in Table 1 summarizes the critical ratios (  ) for an acceptable error = |normalized discrepancy| < 0.01.These   ratios are the numbers highlighted by the arrows in Figures 4 and 5. Column (b) shows the corresponding critical element size ℎ × ℎ to be used for a case with  = 1 m (for isotropic cases,  denotes   or   ; for layered cases,  denotes   ).The use of element size > critical size will lead to |normalized discrepancy| > 0.01.Column (c) summarizes   for an acceptable error = |normalized discrepancy| < 0.05, and Column (d) shows the corresponding critical element size.

Comparison with the Conclusions for Mobilized Shear
Strength.Ching and Phoon [20] 4 and 5 (grey asterisks, reproduced from [20]).For isotropic cases, the absolute normalized discrepancies for Young's modulus are significantly less than those for shear strength (see Figure 4).The large discrepancies for shear strength are due to the fact that the mobilized

Midpoint Strategy
Another possible strategy of discretizing the continuous random field is to adopt a rather naïve midpoint (MP) strategy [16] of assigning E values to elements.In contrast to ELA, the E sample value of the continuous random field at the centroid of each element is adopted.With the MP strategy, the same RFEA with variable element sizes is conducted.However, the process of successive averaging shown in Figure 2 is not adopted; the E value at the centroid of each element is taken, regardless of Sizes 1∼8.Figures 6 and 7 show the normalized discrepancies for isotropic and layered cases, respectively.These two figures should be compared to Figures 4 and 5 for element-level averaging (ELA).
In general, the normalized discrepancies for MP are noticeably larger than those for ELA.This is particularly true for isotropic cases (Figure 4 versus Figure 6).This implies that ELA is a more effective strategy of discretizing the random field than MP, especially for isotropic cases.This observation follows naturally from the fact that  eff for isotropic cases can be well approximated as the geometric average.For comparison, Figures 6 and 7 also show the normalized discrepancies for shear strength (reproduced from [20]).Column (a) in Table 2 summarizes the critical ratios (  ) for an acceptable error = 0.01 for the MP strategy, and Column (b) shows the corresponding critical element size ℎ × ℎ to be used for a case with  = 1 m.Column (c) summarizes   for an acceptable error = 0.05, and Column (d) shows the corresponding critical element size.

Squared Exponential Autocorrelation Model
All the above discussions are for random fields simulated by the single exponential model in (1).In this section, the squared exponential (QExp) model in ( 2) is studied.Table 3 shows the critical ratios (  ) for the ELA strategy, as well as the corresponding critical element sizes.Table 4 shows the  1 and 2), the requirement for QExp is less strict;   for QExp are much smaller than those for SExp.This is because there are significant local oscillations in the SExp random fields, whereas such local oscillations are absent for the QExp random fields [20].

Concluding Remarks
(1) The required mesh size to achieve certain accuracy depends on the discretization strategy (ELA or MP), spatial variability pattern (isotropic or anisotropic), displacement-controlled loading direction (parallel or perpendicular to layers), and autocorrelation model (SExp or QExp).( 2) The element-level averaging (ELA) discretization strategy seems to work quite well for isotropic random fields of Young's modulus.Relatively coarse element can be adopted without inducing large discrepancy between the discretized solution ( eff FEA ) and the continuous solution ( eff ).For instance, if   =   = 1 m, element size smaller than 0.2 m × 0.2 m (  = 5 in Table 1) yields |normalized discrepancy| < 0.01.This is primarily because effective Young's modulus of a domain can be satisfactorily approximated by its geometric average [2,25,26].For layered cases, the required mesh size is slightly smaller (  = 5∼7.8).In general, RFEA with element size < /5 will only induce little discrepancy from the continuous solution if ELA is adopted.In contrast, the midpoint (MP) discretization strategy seems to work relatively poorly for Young's modulus, inducing larger discrepancy to the continuous solution.The above conclusion (ELA outperforms MP) is very different from the conclusion given in [20] for shear strength, because the mobilized shear strength cannot be satisfactorily approximated as the spatial average over a prescribed domain.
(3) The mesh size requirement for the QExp model is less strict than that for the SExp model.A coarser mesh can be used under the QExp model without inducing larger discrepancy to the continuous solution.This conclusion is similar to that for shear strength [20].(4) The observations drawn from this study should apply to cases with simple and uniform stress states.
However, realistic geotechnical problems are mostly subjected to nonuniform stress states, for example, shallow foundation, embankment, slopes, and excavations.Research is ongoing to ascertain the validity of these observations to more realistic and more complex stress states.

Figure 1 :
Figure 1: Spatially variable soil mass under consideration, modeled by FEA mesh.

Figure 2 :
Figure 2: The process of taking averages over 2 × 2 = 4 elements from Size 1 to Size 2.

Figure 4 (
a) shows the normalized discrepancy in sample mean, defined as [(sample mean of   ) - normalized discrepancy     < 0.01 or 0.05
RFEA will be denoted by    and    , where the superscript "RFEA" and subscript "eff" are omitted for brevity.Hence,  ,eff RFEA and  ,eff RFEA for Size 1 mesh (finest mesh) are denoted by   1 and   1 hereon.Two types of spatial variability are considered: (a)

Table 1 :
Critical   ratios and critical element sizes (ELA + SExp)., defined as [(sample STD of   ) -(sample STD of  1 )]/(sample STD of  1 ).The horizontal axis is dimensionless /(element size).The sample mean and sample STD are obtained with    and    samples from 100 random field realizations.Because there are eight mesh sizes, there are eight data points with symbol "◻" in Figure studied the RFEA mesh size effect on the mobilized (overall) shear strength (  ).