Abstract
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.
1. 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–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 two-dimensional 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 () 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 is inevitable because 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 is quantified as follows: [(statistics for discretized ) – (statistics for continuous )]/(statistics for continuous ). 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 can be minimized by adopting very small elements in RFEA, because the difference between 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: what the relationship between the discrepancy and element size is, whether this relationship would change with factors such as the degree of random field anisotropy and loading direction, and 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.
2. 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 (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 and inherent coefficient of variation (COV) , where are the horizontal and vertical coordinates. Therefore, 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 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 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 with mean and variance :where Re[·] denotes the real part of the enclosed complex number; and are independent zero-mean normal random variables with variance given by [15]where and . Besides simulating the point process of the normal random field , the FSM is also able to directly simulate the arithmetic averages of the normal random field 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 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 () because the impact of the spatial variability of Poisson’s ratio is insignificant [2, 25, 26].
2.2. Finite Element Model
The square domain (an elementary soil mass) is modeled by the RFEA mesh shown in Figure 1. The commercial code ABAQUS is adopted. Four-node plane strain elements with full integration (CPE4) are used. Each element follows isotropic elasticity with E = ELA value and with . 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 x-displacement 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 , 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 and (in the global sense, not strain at a point), it can be shown that effective (overall) Young’s modulus for the x direction isIt should be emphasized that is related to the overall 1D compression response of the soil mass in the x direction under . 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 (, ) samples.
3. Random Finite Element Analyses with Variable Element Sizes
Discrepancy between the continuous solution of and any discretized solution will be present. Let the continuous solution be denoted by and let the RFEA solution be denoted by . In this section, detailed analyses will be carried out to understand the magnitude of the normalized discrepancy between and and its relationship with the element size. The RFEA with variable element sizes are taken: eight element sizes corresponding to , and , respectively, denoted by Size 1, Size , and Size 8. Size 1 mesh is the finest with element size = 0.1 × 0.1. There are therefore 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 and . Size 2 mesh uses Size 1 mesh but conducts geometric averaging over elements in Size 1 to obtain the E values for Size 2 (see Figure 2). Size 2 mesh is then subjected to displacement-controlled compression loadings to obtain and . The same procedure is repeated to progressively obtain the meshes for Size 3, Size 4, , and Size 8. For Size mesh, the resulting and will be denoted by and , where the superscript “RFEA” and subscript “eff” are omitted for brevity. Hence, and for Size 1 mesh (finest mesh) are denoted by and hereon. Two types of spatial variability are considered: (a) (isotropic random field) and (b) , (layered random field). Note that with the single exponential model in (1), cases with are not truly isotropic. The term “isotropic” is for notational simplicity.
The normalized discrepancy is defined as [(statistics of ) − (statistics of )]/(statistics of ): it quantifies the discrepancy between the statistics of and the statistics of . Note that serves as the comparison baseline because it is assumed to be close to the continuous solution . 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).
3.1. Isotropic Cases
Figure 4(a) shows the normalized discrepancy in sample mean, defined as [(sample mean of ) (sample mean of )]/(sample mean of ) for , and 100 (isotropic cases). refers to effective Young’s modulus produced by a Size N mesh. Figure 4(b) shows the normalized discrepancy in sample standard deviation (STD), defined as [(sample STD of ) (sample STD of )]/(sample STD of ). 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 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.
(a) Discrepancy in sample mean for
(b) Discrepancy in sample STD for
For isotropic cases, the discrepancy in the sample mean is always positive (i.e., sample mean greater than sample mean). This implies that the use of a coarse mesh in place of a fine mesh will overpredict (unconservative). The discrepancy in the sample STD is also positive (i.e., sample STD greater than 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, 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 , 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 . This ratio of 5 is called the critical ratio () for an acceptable error = 0.01.
3.2. 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 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, 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 , δ/(element size) needs to be greater than 5 ().
(a) Discrepancy in sample mean for
(b) Discrepancy in sample STD for
(c) Discrepancy in sample mean for
(d) Discrepancy in sample STD for
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 (Figures 5(a) and 5(b)). The discrepancy in the sample mean is always positive (i.e., sample mean larger than sample mean). Ching et al. [26] confirmed that for layered cases 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 , δ/(element size) needs to be greater than 7.8 ().
Column (a) in Table 1 summarizes the critical ratios () for an acceptable error . 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 m (for isotropic cases, δ denotes or ; for layered cases, δ denotes ). The use of element size > critical size will lead to . Column (c) summarizes for an acceptable error , and Column (d) shows the corresponding critical element size.
3.3. Comparison with the Conclusions for Mobilized Shear Strength
Ching and Phoon [20] studied the RFEA mesh size effect on the mobilized (overall) shear strength (). They also produced the normalized discrepancy relationships: [(sample mean of ) – (sample mean of )]/(sample mean of ) versus δ/(element size) as well as [(sample STD of ) – (sample STD of )]/(sample STD of ) versus δ/(element size), where denotes the mobilized shear strength simulated by Size N mesh subjected to a compression load. The normalized discrepancies for isotropic and layered cases with are plotted in Figures 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 shear strength cannot be well approximated by any spatial average over a prescribed domain [28, 29]. This is in contrast with effective Young’s modulus () that can be well approximated by the geometric average for isotropic cases. For layered cases, the absolute normalized discrepancies for Young’s modulus and shear strength are comparable (see Figure 5).
Therefore, in order to achieve certain accuracy, for example, , in isotropic cases, the mesh size requirement for Young’s modulus is less strict than that for shear strength. Namely, a coarser mesh can be used for Young’s modulus. However, this conclusion is only limited to isotropic cases, not to layered cases.
4. 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).
(a) Discrepancy in sample mean for
(b) Discrepancy in sample STD for
(a) Discrepancy in sample mean for
(b) Discrepancy in sample STD for
(c) Discrepancy in sample mean for
(d) Discrepancy in sample STD for
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 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 m. Column (c) summarizes for an acceptable error = 0.05, and Column (d) shows the corresponding critical element size.
5. 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 results for the MP strategy. Compared to the SExp results (Tables 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].
6. 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 () and the continuous solution (). For instance, if m, element size smaller than 0.2 m × 0.2 m ( in Table 1) yields . 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This study was under the sponsorship of the Ministry of Science and Technology (MOST) of Taiwan, under Project 103-2221-E-002-129-MY3. The first author would like to acknowledge such a gracious support from MOST. The authors would also like to thank Professor K. K. Phoon for his valuable comments during the writing process.