Finite Element Computations

A recent finite element method for digital images [26] was used to generate data with which to check the results of the D-EMT. Of course, if one is going to check the D-EMT results with the finite element method results, one must ask: how accurate are the finite element method results themselves? Since these results are for concentrated, random systems, there are no exact analytical data against which to check the numerical results. Fortunately, by careful consideration of the possible sources of errors, one can establish the accuracy of the numerical data.

These numerical computations are carried out by first generating the random microstructure desired by building a digital image using square pixels on a 2-D square lattice or cubic pixels on a 3-D cubic lattice. Each pixel is then considered to be a bi-linear element (2-D) or tri-linear element (3-D) [25], so that the entire digital lattice is treated as the finite element mesh. The elastic displacements are linearly extrapolated across the pixels, which is why the the pixels are called bi-linear (2-D) and tri-linear (3-D). The elastic equations are written as a variational principle in the elastic energy, which is then minimized over the digital lattice. The effective moduli are usually defined by a stress average, although they could be defined by an energy average [25].

three main sources of error: (1) finite size effect, (2) digital resolution, and (3) statistical variation [39].

The finite size effect comes about because any given digital image, even with periodic boundary conditions, can only represent a small part of a large random solid. Here we are thinking of inclusions embedded in a matrix. There can be errors induced if the sample is not large enough to possess enough inclusions to be statistically representative. This sampling error can be assessed by running several different size samples, and seeing whether the results change between system sizes. When assessing this source of error, the inclusion size, in terms of pixels, stays the same, so that samples having larger lattices contain more inclusions.

The digital resolution error comes from using square or cubic pixels to represent the inclusions. Even if the inclusions had the same shape as the digital lattice, there would still be a resolution error since one is representing continuum equations with a digital lattice. The size of this error can be checked by holding the number of inclusions constant, and varying the size of the lattice so that there are more or fewer pixels per unit length.

The statistical variation error source simply comes about because the systems under consideration are random ones. For a given concentration of inclusions, there are many ways in which the inclusions might be randomly arranged. Each arrangement will have somewhat different elastic moduli, in general. The size of this error source can be assessed by computing the elastic moduli of several different realizations of the same system (same size lattice, same pixel size and number of inclusions).

In what follows, all the systems shown in the figures were prepared by random sequential adsorption [40], with the largest particles placed first, and then in descending order in diameter. The inner and outer particles were not allowed to overlap with any other inner or outer particle. In the placement process, a trial center was picked at random. If a particle centered at this point did not overlap any other particle, then it was placed and a new site chosen at random. If the particle did overlap another particle, then the site was abandoned and a new site was chosen at random.

Figure 2 shows the 2-D 1000² systems that were used to check the D-EMT results. There were three sizes of inclusions. On a 1000² size digital lattice, these had outer-inner diameters, in pixels, of 121-99, 91-69, and 61-39. Experience with many previous results has shown that having the diameter of the largest inclusion less than one-eighth of the size of the unit cell makes any finite size effects negligible. Holding the number of inclusions fixed to assess the digital resolution error, calculations were also made at sizes of 500² and 2000². There was only 1 % - 2 % variation among the different sizes, so that the size used for all the runs was 1000². The statistical variation for 1000² size systems was very small, and so was neglected.

Figure 2: Gray-scale pictures of the 2-D models used to test the D-EMT predictions (non-overlapping shells). Using concrete terminology, the dark gray is the cement paste matrix, the middle gray is the ITZ regions, and the lightest gray phase is the inclusion or aggregate phase. There are three sizes of aggregates in this picture, and the matrix area fraction is 0.3 (left) and 0.5 (right).

The inner particle had Young's modulus E₃ = 5.0 and Poisson's ratio $\nu_3 = 0.2$₃ = 0.2, the matrix had E₁ = 1.0 and $\nu_1 = 0.3$₁ = 0.3, while the shell had a Young's modulus E₂ that ranged from 0.1 to about 10.0 and a Poisson's ratio of $\nu_2 = 0.3$₂ = 0.3. Two systems were chosen, with matrix area fractions of approximately 0.3 and 0.5. Exact details of the microstructures considered are displayed in Tables 1 and 2. These details are given for future reference, as these highly accurate results for random, non-dilute systems are practically unique.

Table 1: Parameters defining the 2-D microstructures. N = number, L = large circles, M = medium circles, S = small circles, and c is area fraction. Phase 3 is the inner circle, phase 1 is the matrix, and phase 2 is the shell material.
NL	NM	NS	c1	c2	c3
36	28	35	0.303	0.275	0.422
26	20	25	0.499	0.198	0.303

Table 2: Areas of circles and volumes of spheres used in pixels and voxels.
Circle diameter	Sphere diameter	Area (pixel #)	Volume (voxel #)
121		11476
99		7668
91		6488
69		3720
61		2912
39		1184
	28		11536
	16		2176
	14		1472
	8		280

In 3-D, there were two systems considered. The first used only one size of composite sphere, with a ratio of outer to inner diameters of 1.75. The ratio of the unit cell size to the particle outer diameter size was chosen to be about 7, which makes finite size errors negligible. The digital resolution effect was analyzed by running the exact same geometry at sizes of 100³, 200³, and 300³. Figure 3 shows four non-consecutive parallel slices of the 300³ system. A full size range was only run for E₂ = 0.1 and 10.0, which were the limits of the shell stiffness. The phase moduli were the same numerically as in 2-D. We found that at E₂ = 0.1 the resolution error was essentially zero, as all three systems gave almost exactly the same answer, within less than 0.1 %. However, at E₂ = 10.0, there was a 3.7 % drop in bulk modulus and 5.6 % drop in shear modulus between the 100³ and the 200³ systems. There was only a 0.5 % further drop in bulk modulus and 0.8 % in shear modulus when going to the 300³ system, so it was decided that the 200³ system would give adequate accuracy for all the shell moduli. For this size system, the outer particle diameter was 28 pixels wide, while the inner diameter was 16 pixels wide. Statistical errors between configurations were negligible.

Figure 3: Four 300³ slices from the 3-D monosize composite sphere model used to test the D-EMT predictions. Using concrete temrminology, the dark gray is the cement paste matrix, the middle gray is the ITZ regions, and the lightest gray phase is the inclusion or aggregate phase. The matrix volume fraction is 0.668.

The second 3-D system was 200³ in size, and had two size spheres, with outer and inner particle diameters of 28-16 and 14-8. Judging from the data obtained on the single-sphere results, this choice should also have given adequate accuracy, though it was not as carefully checked as were all the other systems. The detailed parameters used for both 3-D systems are presented in Tables 3 and 2. Figure 4 shows four parallel non-consecutive slices of the 200³ system used.

Table 3: Parameters defining the 3-D microstructures, which had either one or two sizes of spheres. The numbers are all for 2003 systems. N = number, L = large spheres, S = small spheres, and c is volume fraction. Phase 3 is the inner sphere, phase 1 is the matrix, and phase 2 is the shell material.
NL	NS	c1	c2	c3
230		0.668	0.269	0.063
220	900	0.517	0.392	0.091

Figure 4: Four 200³ slices from the 3-D two-size sphere model used to test the D-EMT predictions. Using concrete terminology, the dark gray is the cement paste matrix, the middle gray is the ITZ regions, and the lightest gray phase is the inclusion or aggregate phase. The matrix volume fraction is 0.517.