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

Authors: Jianye Ching; Yu-Gang Hu
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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.


---

## Body

## 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 (E
e
f
f) 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 forE
e
f
f is inevitable because E
e
f
f 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 E
e
f
f is quantified as follows: [(statistics for discretized E
e
f
f) – (statistics for continuous E
e
f
f)]/(statistics for continuous E
e
f
f). 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 E
e
f
f can be minimized by adopting very small elements in RFEA, because the difference between E
e
f
f 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.


## 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 L
×
L
=
12.8
×
12.8 (Figure 1). No unit is specified for the dimensionL 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 E
x
,
z, is simulated as a stationary lognormal random field with inherent mean =
μ
=
100 and inherent coefficient of variation (COV) =
V
=
0.5, where (
x
,
z
) are the horizontal and vertical coordinates. Therefore, ln
⁡
E
x
,
z 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
⁡
E
x
,
z between two locations with horizontal distance =
Δ
x and vertical distance =
Δ
z, the single exponential (SExp) autocorrelation model is considered [12, 23]:(1)
ρ
Δ
x
,
Δ
z
=
exp
⁡
-
2
Δ
x
δ
x
-
2
Δ
z
δ
z
,where δ
x and δ
z are the horizontal and vertical SOFs, respectively, for the ln
⁡
E
x
,
z random field. Another popular autocorrelation model is the squared exponential (QExp) model:(2)
ρ
Δ
x
,
Δ
z
=
exp
⁡
-
π
Δ
x
δ
x
2
-
π
Δ
z
δ
z
2
.It is clear that the correlation decreases as ∆
x and ∆
z 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.Figure 1
Spatially variable soil mass under consideration, modeled by FEA mesh.The first author has developed the Fourier series method (FSM) for simulating stationary normal random fields [15, 24]. A 2D stationary lognormal random field E
(
x
,
z
) can be simulated by taking the exponential of a 2D stationary normal random field ln
⁡
E
x
,
z with mean =
λ
=
ln
⁡
[
μ
/
(
1
+
V
2
)
0.5
] and variance =
ξ
2
=
ln
⁡
(
1
+
V
2
):(3)
E
x
,
z
=
exp
⁡
λ
+
R
e
∑
m
=
-
∞
∞
∑
n
=
-
∞
∞
a
m
n
+
i
b
m
n
exp
⁡
i
2
m
π
x
L
x
+
i
2
n
π
z
L
z
,where Re[·] denotes the real part of the enclosed complex number; a
m
n and b
m
n are independent zero-mean normal random variables with variance σ
m
n
2 given by [15](4)
σ
m
n
2
=
ξ
2
q
x
q
z
×
1
-
e
x
p
-
q
x
-
1
m
1
+
m
2
π
2
/
q
x
2
×
1
-
e
x
p
-
q
z
-
1
n
1
+
n
2
π
2
/
q
z
2
S
E
x
p
ξ
2
q
x
q
z
×
e
x
p
-
m
2
π
q
x
2
×
e
x
p
-
n
2
π
q
z
2
Q
E
x
p
,where q
x
=
L
/
δ
x and q
z
=
L
/
δ
z. Besides simulating the point process of the normal random field ln
⁡
[
E
(
x
,
z
)
], the FSM is also able to directly simulate the arithmetic averages of the normal random field ln
⁡
[
E
(
x
,
z
)
] 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
⁡
[
E
(
x
,
z
)
] over that element. By doing this, the local average for E
(
x
,
z
) 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].


### 2.2. Finite Element Model

The12.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. Four-node plane strain elements with full integration (CPE4) are used. Each element follows isotropic elasticity withE = ELA value and with ν
=
0.3. For each realization of theE random field, two RFEAs that simulate an overall 1D compression in thex orz directions are conducted. Let us consider thex 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 uniformx-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 thex direction. The other strain components, ε
z and ε
x
z, 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 thex direction, denoted by E
x
,
e
f
f, 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” σ
x is equal to the arithmetic average σ
x over Side 1. Because the overall strains are ε
x
=
0.1
/
12.8 and ε
z
=
ε
x
z
=
0 (in the global sense, not strain at a point), it can be shown that effective (overall) Young’s modulus for thex direction is(5)
E
x
,
e
f
f
=
1
+
ν
×
1
-
2
ν
1
-
ν
×
overall
σ
x
overall
ε
x
=
1.3
×
0.4
0.7
×
0.1
/
12.8
×
overall
σ
x
=
95.086
×
overall
σ
x
.It should be emphasized that E
x
,
e
f
f is related to the overall 1D compression response of the soil mass in thex direction under ν
=
0.3. A similar displacement-controlled loading is applied in thez direction over the same realization of theE random field. Therefore, eachE random field realization produces a set of (E
x
,
e
f
f, E
z
,
e
f
f) samples.


## 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 L
×
L
=
12.8
×
12.8 (Figure 1). No unit is specified for the dimensionL 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 E
x
,
z, is simulated as a stationary lognormal random field with inherent mean =
μ
=
100 and inherent coefficient of variation (COV) =
V
=
0.5, where (
x
,
z
) are the horizontal and vertical coordinates. Therefore, ln
⁡
E
x
,
z 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
⁡
E
x
,
z between two locations with horizontal distance =
Δ
x and vertical distance =
Δ
z, the single exponential (SExp) autocorrelation model is considered [12, 23]:(1)
ρ
Δ
x
,
Δ
z
=
exp
⁡
-
2
Δ
x
δ
x
-
2
Δ
z
δ
z
,where δ
x and δ
z are the horizontal and vertical SOFs, respectively, for the ln
⁡
E
x
,
z random field. Another popular autocorrelation model is the squared exponential (QExp) model:(2)
ρ
Δ
x
,
Δ
z
=
exp
⁡
-
π
Δ
x
δ
x
2
-
π
Δ
z
δ
z
2
.It is clear that the correlation decreases as ∆
x and ∆
z 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.Figure 1
Spatially variable soil mass under consideration, modeled by FEA mesh.The first author has developed the Fourier series method (FSM) for simulating stationary normal random fields [15, 24]. A 2D stationary lognormal random field E
(
x
,
z
) can be simulated by taking the exponential of a 2D stationary normal random field ln
⁡
E
x
,
z with mean =
λ
=
ln
⁡
[
μ
/
(
1
+
V
2
)
0.5
] and variance =
ξ
2
=
ln
⁡
(
1
+
V
2
):(3)
E
x
,
z
=
exp
⁡
λ
+
R
e
∑
m
=
-
∞
∞
∑
n
=
-
∞
∞
a
m
n
+
i
b
m
n
exp
⁡
i
2
m
π
x
L
x
+
i
2
n
π
z
L
z
,where Re[·] denotes the real part of the enclosed complex number; a
m
n and b
m
n are independent zero-mean normal random variables with variance σ
m
n
2 given by [15](4)
σ
m
n
2
=
ξ
2
q
x
q
z
×
1
-
e
x
p
-
q
x
-
1
m
1
+
m
2
π
2
/
q
x
2
×
1
-
e
x
p
-
q
z
-
1
n
1
+
n
2
π
2
/
q
z
2
S
E
x
p
ξ
2
q
x
q
z
×
e
x
p
-
m
2
π
q
x
2
×
e
x
p
-
n
2
π
q
z
2
Q
E
x
p
,where q
x
=
L
/
δ
x and q
z
=
L
/
δ
z. Besides simulating the point process of the normal random field ln
⁡
[
E
(
x
,
z
)
], the FSM is also able to directly simulate the arithmetic averages of the normal random field ln
⁡
[
E
(
x
,
z
)
] 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
⁡
[
E
(
x
,
z
)
] over that element. By doing this, the local average for E
(
x
,
z
) 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].


## 2.2. Finite Element Model

The12.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. Four-node plane strain elements with full integration (CPE4) are used. Each element follows isotropic elasticity withE = ELA value and with ν
=
0.3. For each realization of theE random field, two RFEAs that simulate an overall 1D compression in thex orz directions are conducted. Let us consider thex 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 uniformx-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 thex direction. The other strain components, ε
z and ε
x
z, 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 thex direction, denoted by E
x
,
e
f
f, 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” σ
x is equal to the arithmetic average σ
x over Side 1. Because the overall strains are ε
x
=
0.1
/
12.8 and ε
z
=
ε
x
z
=
0 (in the global sense, not strain at a point), it can be shown that effective (overall) Young’s modulus for thex direction is(5)
E
x
,
e
f
f
=
1
+
ν
×
1
-
2
ν
1
-
ν
×
overall
σ
x
overall
ε
x
=
1.3
×
0.4
0.7
×
0.1
/
12.8
×
overall
σ
x
=
95.086
×
overall
σ
x
.It should be emphasized that E
x
,
e
f
f is related to the overall 1D compression response of the soil mass in thex direction under ν
=
0.3. A similar displacement-controlled loading is applied in thez direction over the same realization of theE random field. Therefore, eachE random field realization produces a set of (E
x
,
e
f
f, E
z
,
e
f
f) samples.


## 3. Random Finite Element Analyses with Variable Element Sizes

Discrepancy between the continuous solution ofE
e
f
f and any discretized solution will be present. Let the continuous solution be denoted by E
e
f
f and let the RFEA solution be denoted by E
e
f
f
R
F
E
A. In this section, detailed analyses will be carried out to understand the magnitude of the normalized discrepancy between E
e
f
f and E
e
f
f
R
F
E
A 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 inx andz directions to obtain E
x
,
e
f
f
R
F
E
A and E
z
,
e
f
f
R
F
E
A. Size 2 mesh uses Size 1 mesh but conducts geometric averaging over 2
×
2
=
4 elements in Size 1 to obtain theE values for Size 2 (see Figure 2). Size 2 mesh is then subjected to displacement-controlled compression loadings to obtain E
x
,
e
f
f
R
F
E
A and E
z
,
e
f
f
R
F
E
A. The same procedure is repeated to progressively obtain the meshes for Size 3, Size 4, …, and Size 8. For Size N mesh, the resulting E
x
,
e
f
f
R
F
E
A and E
z
,
e
f
f
R
F
E
A will be denoted by E
x
N and E
z
N, where the superscript “RFEA” and subscript “eff” are omitted for brevity. Hence, E
x
,
e
f
f
R
F
E
A and E
z
,
e
f
f
R
F
E
A for Size 1 mesh (finest mesh) are denoted by E
x
1 and E
z
1 hereon. Two types of spatial variability are considered: (a) δ
x
=
δ
z
=
δ (isotropic random field) and (b) δ
x
=
∞,δ
z
=
δ (layered random field). Note that with the single exponential model in (1), cases with δ
x
=
δ
z are not truly isotropic. The term “isotropic” is for notational simplicity.Figure 2
The process of taking averages over2
×
2
=
4 elements from Size 1 to Size 2.The normalized discrepancy is defined as [(statistics ofE
N) − (statistics of E
1)]/(statistics of E
1): it quantifies the discrepancy between the statistics of E
N and the statistics of E
1. Note that E
1 serves as the comparison baseline because it is assumed to be close to the continuous solution E
e
f
f. 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).Figure 3
A logic block diagram for the main task.


### 3.1. Isotropic Cases

Figure4(a) shows the normalized discrepancy in sample mean, defined as [(sample mean of E
N) – (sample mean of E
1)]/(sample mean of E
1) for δ
=
1,5
,
10, and 100 (isotropic cases). E
N refers to effective Young’s modulus produced by a SizeN mesh. Figure 4(b) shows the normalized discrepancy in sample standard deviation (STD), defined as [(sample STD of E
N) – (sample STD of E
1)]/(sample STD of E
1). The horizontal axis is dimensionless δ/(element size). The sample mean and sample STD are obtained with E
x
N and E
z
N 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 E
x
N and E
z
N are the same. Therefore, only the results for E
x
N are presented in Figure 4.Figure 4
Normalized discrepancies for isotropic cases (ELA + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
NFor isotropic cases, the discrepancy in the sample mean is always positive (i.e.,E
x
N sample mean greater than E
x
1 sample mean). This implies that the use of a coarse mesh in place of a fine mesh will overpredict E
e
f
f (unconservative). The discrepancy in the sample STD is also positive (i.e., E
x
N sample STD greater than E
x
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, E
e
f
f 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 E
e
f
f, 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 (r
c) for an acceptable error = 0.01.


### 3.2. Layered Cases

For layered cases, Figures5(a) and 5(b) show the normalized discrepancies in sample mean and STD for E
x
N, effective Young’s modulus when the elementary soil mass is subjected tox direction displacement-controlled loading (parallel to layers). In contrast to isotropic cases, the discrepancy in the sample mean is always negative (i.e., E
x
N sample mean smaller than E
x
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, E
x
,
e
f
f 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 (r
c
=
5).Figure 5
Normalized discrepancies for layered cases (ELA + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
N
(c)
Discrepancy in sample mean forE
z
N
(d)
Discrepancy in sample STD forE
z
NFigures5(c) and 5(d) show the normalized discrepancies in sample mean and STD for E
z
N (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 E
x
,
e
f
f (Figures 5(a) and 5(b)). The discrepancy in the sample mean is always positive (i.e., E
z
N sample mean larger than E
z
1 sample mean). Ching et al. [26] confirmed that for layered cases E
z
,
e
f
f 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 (r
c
=
7.8).Column (a) in Table1 summarizes the critical ratios (r
c) for an acceptable error =
|
normalized discrepancy
|
<
0.01. These r
c ratios are the numbers highlighted by the arrows in Figures 4 and 5. Column (b) shows the corresponding critical element size h
×
h to be used for a case with δ
=
1 m (for isotropic cases, δ denotes δ
x or δ
z; for layered cases, δ denotes δ
z). The use of element size > critical size will lead to |
normalized discrepancy
|
>
0.01. Column (c) summarizes r
c for an acceptable error =
|
normalized discrepancy
|
<
0.05, and Column (d) shows the corresponding critical element size.Table 1
Criticalr
c ratios and critical element sizes (ELA + SExp).
Case
⁢0.01-tolerance
⁢0.05-tolerance
r
c (a)
Critical element size
h
×
h for δ = 1 m(b)
r
c (c)
Critical element size
h
×
h for δ = 1 m(d)
Isotropic case
5
0.2 m
0.078
12.82 m
Layered caseE
x
,
e
f
f
5
0.2 m
1
1 m
Layered caseE
z
,
e
f
f
7.8
0.13 m
2
0.5 m


### 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 (τ
f). They also produced the normalized discrepancy relationships: [(sample mean of τ
f
N) – (sample mean of τ
f
1)]/(sample mean of τ
f
1) versus δ/(element size) as well as [(sample STD of τ
f
N) – (sample STD of τ
f
1)]/(sample STD of τ
f
1) versus δ/(element size), where τ
f
N denotes the mobilized shear strength simulated by SizeN mesh subjected to a compression load. The normalized discrepancies for isotropic and layered cases with δ
=
10 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 (E
e
f
f) 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,|
normalized discrepancy
|
<
0.01, 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.


## 3.1. Isotropic Cases

Figure4(a) shows the normalized discrepancy in sample mean, defined as [(sample mean of E
N) – (sample mean of E
1)]/(sample mean of E
1) for δ
=
1,5
,
10, and 100 (isotropic cases). E
N refers to effective Young’s modulus produced by a SizeN mesh. Figure 4(b) shows the normalized discrepancy in sample standard deviation (STD), defined as [(sample STD of E
N) – (sample STD of E
1)]/(sample STD of E
1). The horizontal axis is dimensionless δ/(element size). The sample mean and sample STD are obtained with E
x
N and E
z
N 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 E
x
N and E
z
N are the same. Therefore, only the results for E
x
N are presented in Figure 4.Figure 4
Normalized discrepancies for isotropic cases (ELA + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
NFor isotropic cases, the discrepancy in the sample mean is always positive (i.e.,E
x
N sample mean greater than E
x
1 sample mean). This implies that the use of a coarse mesh in place of a fine mesh will overpredict E
e
f
f (unconservative). The discrepancy in the sample STD is also positive (i.e., E
x
N sample STD greater than E
x
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, E
e
f
f 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 E
e
f
f, 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 (r
c) for an acceptable error = 0.01.


## 3.2. Layered Cases

For layered cases, Figures5(a) and 5(b) show the normalized discrepancies in sample mean and STD for E
x
N, effective Young’s modulus when the elementary soil mass is subjected tox direction displacement-controlled loading (parallel to layers). In contrast to isotropic cases, the discrepancy in the sample mean is always negative (i.e., E
x
N sample mean smaller than E
x
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, E
x
,
e
f
f 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 (r
c
=
5).Figure 5
Normalized discrepancies for layered cases (ELA + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
N
(c)
Discrepancy in sample mean forE
z
N
(d)
Discrepancy in sample STD forE
z
NFigures5(c) and 5(d) show the normalized discrepancies in sample mean and STD for E
z
N (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 E
x
,
e
f
f (Figures 5(a) and 5(b)). The discrepancy in the sample mean is always positive (i.e., E
z
N sample mean larger than E
z
1 sample mean). Ching et al. [26] confirmed that for layered cases E
z
,
e
f
f 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 (r
c
=
7.8).Column (a) in Table1 summarizes the critical ratios (r
c) for an acceptable error =
|
normalized discrepancy
|
<
0.01. These r
c ratios are the numbers highlighted by the arrows in Figures 4 and 5. Column (b) shows the corresponding critical element size h
×
h to be used for a case with δ
=
1 m (for isotropic cases, δ denotes δ
x or δ
z; for layered cases, δ denotes δ
z). The use of element size > critical size will lead to |
normalized discrepancy
|
>
0.01. Column (c) summarizes r
c for an acceptable error =
|
normalized discrepancy
|
<
0.05, and Column (d) shows the corresponding critical element size.Table 1
Criticalr
c ratios and critical element sizes (ELA + SExp).
Case
⁢0.01-tolerance
⁢0.05-tolerance
r
c (a)
Critical element size
h
×
h for δ = 1 m(b)
r
c (c)
Critical element size
h
×
h for δ = 1 m(d)
Isotropic case
5
0.2 m
0.078
12.82 m
Layered caseE
x
,
e
f
f
5
0.2 m
1
1 m
Layered caseE
z
,
e
f
f
7.8
0.13 m
2
0.5 m


## 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 (τ
f). They also produced the normalized discrepancy relationships: [(sample mean of τ
f
N) – (sample mean of τ
f
1)]/(sample mean of τ
f
1) versus δ/(element size) as well as [(sample STD of τ
f
N) – (sample STD of τ
f
1)]/(sample STD of τ
f
1) versus δ/(element size), where τ
f
N denotes the mobilized shear strength simulated by SizeN mesh subjected to a compression load. The normalized discrepancies for isotropic and layered cases with δ
=
10 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 (E
e
f
f) 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,|
normalized discrepancy
|
<
0.01, 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 assigningE values to elements. In contrast to ELA, theE 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; theE 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).Figure 6
Normalized discrepancies for isotropic cases (MP + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
NFigure 7
Normalized discrepancies for layered cases (MP + SExp).
(a)
Discrepancy in sample mean forE
x
N
(b)
Discrepancy in sample STD forE
x
N
(c)
Discrepancy in sample mean forE
z
N
(d)
Discrepancy in sample STD forE
z
NIn general, the normalized discrepancies for MP are noticeably larger than those for ELA. This is particularly true for isotropic cases (Figure4 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 E
e
f
f 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 (r
c) for an acceptable error = 0.01 for the MP strategy, and Column (b) shows the corresponding critical element size h
×
h to be used for a case with δ
=
1 m. Column (c) summarizes r
c for an acceptable error = 0.05, and Column (d) shows the corresponding critical element size.Table 2
Criticalr
c ratios and critical element sizes (MP + SExp).
Case
⁢0.01-tolerance
⁢0.05-tolerance
r
c (a)
Critical element size
h
×
h for δ = 1 m(b)
r
c (c)
Critical element size
h
×
h for δ = 1 m(d)
Isotropic case
50
0.02 m
3.125
0.32 m
Layered caseE
x
,
e
f
f
20
0.05 m
2.5
0.4 m
Layered caseE
z
,
e
f
f
20
0.05 m
2.5
0.4 m


## 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 (r
c) 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; r
c 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].Table 3
Criticalr
c ratios and critical element sizes (ELA + QExp).
Case
⁢0.01-tolerance
⁢0.05-tolerance
r
c (a)
Critical element size
h
×
h for δ = 1 m(b)
r
c (c)
Critical element size
h
×
h for δ = 1 m(d)
Isotropic case
2.5
0.4 m
0.078
12.82 m
Layered caseE
x
,
e
f
f
2.5
0.4 m
0.78
1.82 m
Layered caseE
z
,
e
f
f
2.5
0.4 m
1.25
0.8 mTable 4
Criticalr
c ratios and critical element sizes (MP + QExp).
Case
⁢0.01-tolerance
⁢0.05-tolerance
r
c (a)
Critical element size
h
×
h for δ = 1 m(b)
r
c (c)
Critical element size
h
×
h for δ = 1 m(d)
Isotropic case
5
0.2 m
1.56
0.64 m
Layered caseE
x
,
e
f
f
2.5
0.4 m
1.25
0.8 m
Layered caseE
z
,
e
f
f
2.5
0.4 m
1.56
0.64 m


## 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 (E
e
f
f
F
E
A) and the continuous solution (E
e
f
f). For instance, if δ
x
=
δ
z
=
1 m, element size smaller than 0.2 m × 0.2 m (r
c
=
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 (r
c = 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.



---

*Source: 8756271-2016-08-03.xml*