Mechanical properties of sintered meso-porous silicon: a numerical model

Because of its optical and electrical properties, large surfaces, and compatibility with standard silicon processes, porous silicon is a very interesting material in photovoltaic and microelectromechanical systems technology. In some applications, porous silicon is annealed at high temperature and, consequently, the cylindrical pores that are generated by anodization or stain etching reorganize into randomly distributed closed sphere-like pores. Although the design of devices which involve this material needs an accurate evaluation of its mechanical properties, only few researchers have studied the mechanical properties of porous silicon, and no data are nowadays available on the mechanical properties of sintered porous silicon. In this work we propose a finite element model to estimate the mechanical properties of sintered meso-porous silicon. The model has been employed to study the dependence of the Young’s modulus and the shear modulus (upper and lower bounds) on the porosity for porosities between 0% to 40%. Interpolation functions for the Young’s modulus and shear modulus have been obtained, and the results show good agreement with the data reported for other porous media. A Monte Carlo simulation has also been employed to study the effect of the actual microstructure on the mechanical properties.


Background
Porous silicon (PSi) has been extensively employed in microelectromechanical systems (MEMS) technology and it has been proposed for some applications in photovoltaics (PV) technology. In MEMS technology, processes to manufacture suspended structures often employ openporosity PSi as a sacrificial layer while in PV technology, stacked layers of sintered PSi with different porosity have been proposed both as buried Bragg reflectors and for layer transfer techniques for the fabrication of thin silicon solar cells [1]. In these applications, an accurate evaluation of PSi mechanical properties is paramount for the device fabrication and performance.
Since its discovery, PSi has been investigated mostly for its optical and electrical properties. Only few researchers investigated the mechanical properties of this material. Characterizations of PSi mechanical properties were performed employing very different techniques, e.g., nanoindentation [2], Brillouin scattering [3], phase velocity scanning [4], and microechography [5], but each of these works dealt with open-porosity PSi.
For some applications such as the layer transfer technique and the buried Bragg reflectors, as-anodized PSi undergoes a thermal treatment which transforms the structure of the silicon layer. When annealed at temperatures higher than 1,000°C, PSi experiences a structural reorganization driven by surface energy minimization which transforms the columnar voids formed by the anodization into closed sphere-like pores. The pores interact and give rise to more complex structures (see Figure 1). To the knowledge of the authors, no studies have been performed about the mechanical properties of the sintered meso-porous silicon. Magoariec and Danescu [6] proposed a numerical model for evaluating the mechanical properties of sintered nano-PSi which they compared with experimental data from literature, but they did not take into account the effect of the random distribution of voids. Li et al. [7] showed that the mechanical properties of porous media depend on the connectivity of the pores http://www.nanoscalereslett.com/content/7/1/597 and, thus, the actual distribution of voids is expected to play a role in the mechanical properties of porous silicon.
In this work we present a finite element model for the evaluation of the overall mechanical properties of sintered PSi. This numerical model includes a random distribution of voids inside a representative volume element (RVE), and it has been exploited to evaluate the upper and lower bounds of the Young's modulus and the shear modulus of PSi as a function of porosity. A Monte Carlo simulation is also presented in this work to evaluate the fluctuation of the upper bounds due to the statistical variation of the microstructure.

Methods
An example of the numerical model employed for the analysis of the mechanical properties of PSi is depicted in Figure 2. The RVE is generated by subtracting a random distribution of spherical voids to a silicon cube. Without loss of generality, we assume that the edges of the cube are oriented along the < 100 > directions of the silicon crystal lattice. Referring to this orientation, the Young's modulus, the Poisson ratio and the shear modulus of silicon employed in the analysis are the ones reported by Masolin et al. [8], i.e. E Si = 130 GPa, ν Si = 0.278 and G Si = 79.6 GPa. The position of the center of the voids and their radii follow a rectangular and a Gaussian distribution, respectively. Since the single void can be centered in any location in the silicon cube, there is a nonnull probability that two or more voids overlap. The RVEs are then discretized by means of simplicial meshes.
To define the overall mechanical properties of PSi, the standard homogenization theory has been employed. In this framework, three different kinds of boundary conditions are usually applied: uniform displacements, uniform tractions, and periodic boundary conditions. If the considered RVE is sufficiently large, these boundary conditions lead to the same overall mechanical response, but such RVE is usually excessively large to be solved by numerical simulations. For smaller RVEs, uniform tractions and uniform displacements boundary conditions are known to evaluate respectively the upper and lower bounds of the stiffness matrix of the RVE, while periodic boundary conditions give a more reasonable estimation of the homogenized stiffness matrix. In this work we want to focus only on the values that bound the mechanical properties and, thus, periodic boundary conditions will not be taken into account.
In the case of uniform displacement boundary conditions, displacements u i have been imposed at the boundaries (∂V ) of the RVE (V ) by the relation: where E ij are the components of the second order macroscopic strain tensor while, in case of uniform tractions boundary conditions, uniform surface loads have been applied at the boundary of the RVE. In the latter case, the components of the macroscopic strain tensor are computed from the average of the displacement field at the boundaries. Since the aspect ratio of the structures inside the RVE is limited and since, at room temperature, silicon behaves as a linear elastic material, nonlinear effects can be neglected until the stress field locally reaches the strength of silicon.
The displacements and the loads imposed at the boundary induce a stress field inside the matrix that can be described by the second order microscopic stress tensor σ lm (x). The second order macroscopic stress tensor lm is evaluated as the average of the microscopic stress field in the RVE: where |V | denotes the volume of the RVE. Once the components of the macroscopic stress tensor are computed, the components of the fourth order stiffness tensor (S ijlm ) can be extracted by the appropriate ratio between E ij and lm , following the relation: Since bulk silicon shows an orthotropic behavior with equivalent directions along < 100 >, a similar behavior is assumed for the RVE. By making this assumption, the equivalent Young's modulus and shear modulus in the < 100 > directions can be computed. It has to be noted that the properties obtained by assuming an orthotropic material automatically degenerate in isotropic condition whenever the RVE is not anisotropic. To study the effect of the actual microstructure on the overall mechanical properties, a Monte Carlo simulation has been performed by generating different realizations with the same statistical distribution of pores positions and radii.

Results and discussion
Using the method presented in the previous section, the upper and the lower bounds of both the Young's modulus and the shear modulus of 500 × 500 × 500 nm 3 PSi cubes is the mechanical parameter as function of porosity , A Si is the mechanical parameter of the matrix, i.e. silicon, and k is the only fitting parameter that has to be tuned in the interpolation. This family of functions is commonly employed for the interpolation of mechanical properties of porous solids [9] and, in this work, they have been employed both for the upper bounds (UBs) and the lower bounds (LBs). The interpolating functions and the relative R 2 values obtained for the Young's modulus are as follows: while for the shear modulus, As the R 2 values suggest, these functions fit well the mechanical properties obtained by the simulations.
The computed values are compared with the Voigt bound that is known to define the theoretical upper bound for the elastic moduli, while the Reuss lower bound is neglected since it is trivially null. As expected from the homogenization theory, the values gathered from simulations are between the theoretical bounds.
In Figure 5, the degree of anisotropy computed as G( ) 2(1+ν( )) E( ) , i.e., the ratio between the computed shear modulus and the one evaluated by accounting silicon as an isotropic material is reported for both the upper bound and the lower bound. Even though the degree of anisotropy tends to the unity and, therefore, to the isotropic behavior, sintered PSi still behaves as an orthotropic material also for relatively high porosities.
Since each realization of the Monte Carlo simulation is generated by keeping constant the statistical distribution of the pores radii and positions, the porosity cannot be also fixed. An a-posteriori analysis reveals that the porosity ranges between 15% and 20% with a mean value of 17.66%. Figure 6 shows the statistics of the upper bound of both the Young's modulus and the shear modulus. Within the legend, the mean values and the variances are reported. To discern the effect of the actual pores distribution from the random change of porosity, the correlation coefficients (r E and r G ) between the logarithm of the two moduli and the logarithm of (1 − ) have been computed. Both the correlation coefficients (r E = 0.9485 and r G = 0.9773) report a strong correlation between the porosity and the moduli. This suggests that, once the statistical pores distribution is fixed, the actual microstructure has secondary effect compared to small variations of porosity.

Conclusions
In this work we presented a finite element model to evaluate the overall mechanical properties of sintered meso-PSi. This model has been employed to characterize the upper and lower bound of the Young's modulus and shear modulus for meso-PSi with porosity between 0% and 40%. The values defined by the simulations can be fitted by interpolation functions that are commonly employed for porous media. This analysis reduces the theoretical bounds on the mechanical properties defined by the Voigt and Reuss limits and provides an indication on the possible Young's modulus and shear modulus of sintered PSi as function of porosity. The values of Young's modulus and shear modulus obtained by simulations are well fitted by interpolation functions that have been already employed for other porous media. This suggests that the model could represent well the actual properties of sintered PSi.
Monte Carlo simulations have also been employed to analyze the effect of the actual microstructure on the upper bound for porosities between 15% and 20%. The results show that the large spread on the values of the Young's modulus and shear modulus is mainly due to the variation of the porosity instead of the variation of the actual microstructure itself.
The obtained values can be employed for the optimization of structures which involve sintered PSi, and the model can be exploited to study sintered PSi with different pores distributions.