Finite Element Method

The finite element method uses a variational formulation of the linear elastic equations, and finds the solution by minimizing the elastic energy via a fast conjugate gradient method. The FEM we use has been especially adapted for periodic rectangular parallelepiped digital images (although they can be used on non-periodic images), and is for linear elasticity only. Each pixel, in 3-D, is taken to be a tri-linear finite element [21]. The digital image is assumed to have periodic boundary conditions. A strain is applied, with the average stress or the average elastic energy giving the effective elastic moduli [18,22]. Details of the theory and copies of the actual programs used are reported in the papers of Garboczi & Day [14] and Garboczi [23].

Given a digital microstructure, the finite element method provides a numerical solution of the elasticity equations. The accuracy is only limited by the finite number of pixels which can be used (around 10⁶ in this study). Preliminary studies indicated that about 100 cells are necessary to properly simulate the macroscopic properties of a cellular solid (which may have many thousands of cells). We generally calculated the properties of five samples at each density and report an average value. The statistical uncertainty in the results is estimated to be less than 10 %. Note that if a foam is regular and periodic, just one measurement on a unit cell is sufficient.

A potentially greater source of error occurs in the finite element method when there are insufficient pixels in a solid region to correctly model continuum elasticity. A useful method of estimating these discretization errors is to compute the properties of regular periodic foams, since in these models there are no sources of statistical error, and there are exact solutions to which to compare the numerical results. The foams we consider have cubic symmetry, which means that the direction dependent elastic properties can be characterized by three independent constants, C₁₁, C ₁₂ and C₄₄ , of the Hooke's law stress-strain tensor [24]. For loading along the x-axis (ie. one of the axes of symmetry) the Young's modulus and Poisson's ratio are

$$\frac{(C_{11}-C_{12})(C_{11}+2C_{12})}{(C_{11}+C_{12})}$$ E₁₀₀ = (7)

$$\nu_{12} \frac{C_{12}}{(C_{11}+C_{12})}.$$ = (8)

The bulk modulus is actually independent of direction and given by K = E₁₀₀ / 3(1 - 2₁₂ ), and the anisotropic shear modulus (for shearing parallel to a symmetry plane) is just C₄₄ . The finite element codes evaluate the C_ijdirectly, but for simplicity we report the engineering constants.

To check the effect of resolution for the finite element method we measured the Young's modulus (E₁₀₀ ) of two tetrakaidecahedral models with edges of thickness 4 pixels and 8 pixels, respectively. We found virtually no difference (< 1 %) in E₁₀₀ indicating that the discretization errors are quite small. However, the absolute value is 15 % greater than that estimated by Simone and Gibson [5] using a specially tailored finite-element grid. Since we have found our FEM to be accurate for many other test cases, the origin of the discrepancy is unclear. In related studies, we have found discretization errors of around 10 %, and we assume this will be true for the random models studied here.