17324
J. Phys. Chem. C 2010, 114, 17324–17332
Simulations of Cracks Supported by Experiments: The Influence of the Film Height and Isotropy on the Geometry of Crack Patterns Johannes Richardi, Anh Tu Ngo, and Marie Paule Pileni* Laboratoire des Mate´riaux Me´soscopiques et Nanome´triques (LM2N), UMR CNRS 7070, UniVersite´ Pierre et Marie Curie, baˆt F, BP 52, 4 Place Jussieu, 75252 Paris Cedex 05, France ReceiVed: May 12, 2010; ReVised Manuscript ReceiVed: August 18, 2010
Here, we show that cracks give valuable information on the height and the isotropy of nanocrystal films. Two-dimensional crack patterns are systematically studied by simulations and experiments varying the height and the anisotropy of the applied stress. The simulations are carried out using a bundle-spring network model, which allows studying the influence of height on crack patterns. For the experiments, a model system made of magnetic γ-Fe2O3 nanocrystals is used, which enables a change in the anisotropy of the stress by applying a magnetic field during the drying process. The crack pattern morphology is investigated by simulations using square, hexagonal, and isotropic spring arrangements and applying an isotropic or unilateral stress. The average crack distance as a function of the film height studied by simulations follows a universal scaling law, confirming the experimental data. We show that this implies that the morphology of the crack pattern does not change with the height. The frequency of crack fragments as a function of the number of their sides or neighbors is analyzed for isotropic crack patterns obtained from the simulations using various spring arrangements. Finally, the theoretical results are compared to the experiments. The observed linear scaling of crack distances with height is in good agreement with the theory. The frequency of crack fragments as a function of the number of sides and neighbors does not significantly change with height, except for very thin layers. The simulations show that the experimental frequencies indicate the isotropy of the nanocrystal films. 1. Introduction Shrinkage-crack patterns arise from the stress that builds up when the contraction of a film due to cooling or drying is restricted by the adhesion to a substrate.1,2 Usually, the crack development is sufficiently slow so that the layer is homogeneous over its thickness and the resulting crack patterns are effectively two-dimensional. Such crack patterns are common in both natural and man-made systems,1,3-5 including dried mud layers, paint or other coatings, glazes on ceramic mugs, water-cornstarch mixtures, etc. They are encountered in a great variety of sizes. Thus, mud cracks are observed on drying puddles of a few centimeters in diameter but also on playas with giant polygons reaching sizes of several hundred meters. The marked similarities of crack patterns in various systems and on a large scale of size indicate that they do not depend on the material and the size. Thus, usually experimentalists observe polygonal patterns consisting mostly of straight cracks meeting at 90°. Only in very thin films, where heterogeneities start to play a role, is there a large fraction of junction angles greater than 90° and wavy cracks appear.3-5 The analysis of crack formation is important for the understanding of the mechanism of drying,6 but the cracks also give valuable information on the structure of films.7 Much effort was done to avoid cracks, for example, by decreasing the evaporation rate.8,9 Besides drying, cracks may also be induced by the degradation of organic materials,10 the application of a high voltage,11 or by the mismatch of thermal expansion coefficients.12 The morphology of crack patterns changes by directional drying.4,13,14 In this way, ladder-like patterns4,13 or helices14 can be fabricated. The morphology of cracks can also be controlled * To whom correspondence
[email protected].
should
be
addressed.
E-mail:
by the boundary conditions during the deposition.15 Surprisingly, the morphology seemed to influence the scaling of crack patterns. Thus, experiment13 and theory16 predicted that the average distance between directional cracks varies nonlinearly with the height, in contrast to the linear relation found in the isotropic case.3,4,17 The formation of 2D crack patterns has been widely studied using spring network models and taking the substrate into account. Let us first describe some of the theoretical work done using this model. The spring network approach was pioneered by, for example, P. Meakin18,19 to explain the cracks observed by drying monolayers of close-packed polystyrene beads. Andersen et al.20-22 proposed a new spring network model, which contains only one parameter, the ratio between the thresholds for block slipping and spring breaking. This model was the basis of the work of Leung and Ne´da,17 who proposed to replace the single spring by a bundle of h springs, where h represents the height of the film. In good agreement with a number of experiments, the simulations by P. Meakin18,19 predict that the first cracks are straight, whereas later crack generations are wavy as the cracking process releases the strain. Andersen et al.20-22 observed that the growth of domains of positive and negative components of the stress field is crucial for the formation of the crack patterns. The linear scaling of the crack distance with height was theoretically shown for isotropic crack patterns.17 In several papers, the influence of local disorder on the formation of the size of crack fragments and the propagation of cracks using spring-network models was studied.23,24 They showed that, for small disorder, the film fragments through crack propagation, whereas in highly disordered systems, the cracks are formed by the coalescence of initially independent point defects. Spring network models were also used to reproduce the properties of
10.1021/jp104329g 2010 American Chemical Society Published on Web 09/28/2010
Simulations of Cracks Supported by Experiments directional cracks.25-27 Thus, they were applied to describe the scaling of cracks in thin nonhomogeneous slabs of clay supported by a deformable substrate.25,26 Other examples are the reproduction of helical-shaped cracks, mentioned above,7,27 and the instability of directional crack propagation in thin glass plates.28-32 It should be noted that scalar electrical models to describe the drying-induced crack formation in analogy to the substratesupported spring network model have been developed.33-36 The evolution of crack patterns with a successive hierarchy of several generations of cracks was also widely studied by other approaches.1,37-39 In particular, Bonn et al.37 have recently developed a theory, which explains that fragments bounded by cracks have, on average, four sides and six neighbors. The propagation of the crack tip and the roughness of the fracture surfaces have also attracted much interest.40-47 The formation of cracks has also been studied by atomistic simulations.48 Recently, we have developed a combined approach using experiments and simulations.49,50 In the experiments, we have proposed a new model system using drying films of maghemite nanocrystals.51 Although cracks in nanocrystal and nanoporous systems have been widely observed,52-54 they were rarely used directly to study crack patterns. With respect to other model systems, such as coffee-water mixtures, the new system has two major advantages. First, the small size of the nanocrystals makes possible the study of crack patterns from micrometer to millimeter scales. Second, the magnetization of the particles enables the fabrication of directional crack patterns by the application of a magnetic field, which is an alternative approach to the directional drying method. In the theoretical part of our previous papers, we have applied the bundle-spring network model17 to study directional and other anisotropic crack patterns.49,50 This implied the use of an about 25-fold larger network with respect to the literature to achieve the same statistical accuracy as in the isotropic case. We have shown both by simulations and by experiments that the average crack distance linearly scales with the height for isotropic and directional patterns.49,50 To correctly analyze the experimental data, it is important to differentiate between primary and secondary crack generations; the latter disappear at smaller film heights. Please note that the crack formation in films of magnetic nanocrystals has also very recently been studied by Pauchard et al.55 The deviation of the crack propagation observed experimentally due to a field applied perpendicular to the crack direction is in good agreement with an analytical theory. This approach allows obtaining the Young’s modulus of the film before cracking. Here, the influence of the film thickness on the detailed geometry of isotropic and directional crack patterns is systematically studied by simulations using three different arrangements of springs. Although previous publications49,50 focused on the average crack distance, we will also investigate the angles between crack junctions and numbers of edges and neighbors of crack fragments. The theoretical predictions will be compared to experimental results. The paper is organized as follows. In section 2, the simulation method will be presented. In particular, the technique to produce new isotropic arrangements of springs will be explained. In section 3.1, the crack patterns, obtained at various heights and using different spring arrangements, will be shown and discussed. The cases of an isotropic and a unilateral stress will be treated. In section 3.2, the scaling of the average crack distance for both stress cases, obtained by simulations using various spring arrangements, will be studied. A new proof of the linear scaling law of crack distances with
J. Phys. Chem. C, Vol. 114, No. 41, 2010 17325 height found in the simulations will be described and discussed. In section 3.3, the geometry of the simulated crack patterns obtained under isotropic stress will be quantitatively analyzed by the numbers of edges and neighbors of the crack fragments. In section 3.4, the quality of the theoretical predictions will be tested by comparison with experiments. 2. Simulation Method The simulation method used to study the variation of crack patterns with the film thickness is based on a recently introduced bundle-spring network model proposed by Leung and Ne´da.17 Within this approach, the layer is represented on a mesoscopic scale by a 2D array of blocks where each neighbor is connected by coil springs. The simulations are carried out for a large array of blocks connected by bundles of springs. This model introduced by Leung and Neda applies to a film the thickness of which is much smaller than the average distance of cracks (in our experiments, the ratio between crack distance and thickness is 5.4). In this case, no cracks are expected parallel to the substrate, which is actually experimentally observed. Starting from a 3D model, we can, in a first approximation, replace the springs perpendicular to the substrate by inelastic links. This leads us directly to our 2D model with bundles of springs. Obviously, the contraction perpendicular to the substrate is ignored. Therefore, this model cannot correctly predict the detailed formation of the cracks, which is expected to go from the top of the film to its bottom. However, the geometry of the cracks perpendicular to the film surface is expected to be correct. This approach is experimentally justified because the cracks go straight from the top to the bottom of the films, as observed in the experimental images in our previous papers.49,50 Two blocks i and j are separated by an initial distance aij g l, where l is the relaxed spring length. The blocks are slightly randomly displaced by Dr b, where |Dr b| < 0.001. These random displacements make producing different initial conditions possible. A change of these random displacements influences the exact positions of the cracks, but not the average properties, such as the average surface area limited by the cracks. The force at a block i is given by a harmonic approximation using Hooke’s law17
b Fi )
b r
∑ hijbFij with bFij ) k(rij - l) rijij
(1)
j
where j denotes the neighboring blocks at a distance rij, k is the spring constant fixed at k ) 1, and the spring length l ) 0.9 is chosen. The dimensionless number hij corresponds to the number of springs between the blocks i and j. Before the simulation starts, every hij is set to the value of h. The harmonic approximation is used here because, in realistic situations, the stress is primarily tensile and the strain is low. Three different arrangements of blocks are used for the simulations: (i) a hexagonal grid of points, as previously used in refs 17, 49, and 50, (ii) a square array where only the closest neighbors are connected by springs, as already used before in refs 20 and 49, and (iii) new isotropic arrangements, the generation of which will be explained later on. In contrast to the third arrangement, the first two ones are characterized by a marked anisotropy due to the periodic alignment of grid points. Because the experiments are usually performed on disordered isotropic materials, the isotropic arrangement is expected to be the best one to use. We decided to try a new arrangement because, as we will show in this article,
17326
J. Phys. Chem. C, Vol. 114, No. 41, 2010
Richardi et al.
the other arrangements were unable to reproduce the exact experimental geometries, for example, the number of sides and neighbors of crack fragments. The generation of isotropic arrangements is explained in detail in the Appendix. Please note that an isotropic or unilateral stress will be applied to all three arrangements in this article. To avoid misunderstanding with the term isotropic used for both the stress and the arrangement, we will always use the combinations “isotropic arrangement, grid, or array” and “isotropic stress or pattern”. At the start of each simulation, a strain is applied to the relaxed array in only one or both directions by choosing side lengths Lx and Ly, larger than the values Lx0 and Ly0 of the relaxed block array, to generate the stress. The initial extension of the arrays is defined by
sx ) (Lx - L0x )/Lx and sy ) (Ly - L0y )/Ly For sx ) sy, isotropic 2D crack patterns are obtained. A unilateral stress (sx > 0 and sy ) 0) produces cracks in the y direction. In the simulation, a block is only moved to a new mechanical equilibrium position when the force Fi at the block is larger than a slipping threshold denoted by Fs. Whenever the force Fij on a spring exceeds a cracking threshold denoted by Fc, a spring is broken (hij f hij - 1). During the simulation, the two thresholds are systematically decreased to induce slipping and cracking. Because both thresholds decrease by the same physical means, the ratio κ ) Fc/Fs remains fixed. The blocks in our model represent mesoscopic elements of the film. Therefore, Fc or Fs are the cohesive forces between parts of the film or between the film and the substrate. To ensure that the results do not depend on the chosen κ value, simulations are carried out for κ ranging from 0.25 to 3.0. The details of the simulation run are explained in the Appendix. Please note that, as explained in the Appendix, for the exact simulation of crack formation under unilateral stress, the number of blocks had to be increased from 10 000 (as used in the literature17) to 250 000 blocks.
Figure 1. SEM (a, b) and photographic (c, d) images of isotropic 2D and directional (1D-2D) crack patterns made of 10 nm γ-Fe2O3 nanocrystals produced by applying, or not applying, a magnetic field (0.2 T) during the evaporation process. (a) Micrometric-sized isotropic crack patterns at the end of the drying process. The average layer height is 2.8 ( 0.02 µm. The solution volume and concentration are 0.2 mL and 2% by weight before drying. (b) Micrometric-sized directional crack patterns obtained at the end of the drying process. The average layer height is 7 ( 0.2 µm. The solution volume and concentration are 0.2 mL and 0.35% by weight before drying. (c) Centimeter-sized isotropic crack patterns with an average height of 1000 ( 20 µm. The solution volume and concentration are 27 mL and 43% by weight before drying. (d) Centimeter-sized directional crack patterns with an average height of 415 ( 5 µm. The solution volume and concentration are 5 mL and 28% by weight before drying.
3. Results and Discussions Figure 1 shows typical crack patterns that form in the experimental systems under isotropic (Figure 1a,c) and unilateral stresses (Figure 1b,d). The unilateral stress is produced due to the application of a magnetic field parallel to the substrate. By varying the initial concentration and the volume of the nanocrystal solution, the film height can be changed from micrometer (Figure 1a,b) to centimeter scales (Figure 1c,d). Please note that, with increasing film thickness (Figure 1d), cracks perpendicular to the direction of the applied field are observed. This indicates that the stress is anisotropic, but not strictly unilateral. The aim of the following simulations is to reproduce quantitatively the geometry of the experimental crack patterns in Figure 1. 3.1. Crack Patterns from Simulations under Isotropic and Unilateral Stresses. Here, we compare the crack patterns obtained by simulations using isotropic, hexagonal, and square arrangements of blocks at various film heights. Figure 2 shows the crack patterns observed at the end of the simulations, when initially an isotropic stress is applied to the block array (sx ) 0.1 and sy ) 0.1). The crack patterns obtained using a hexagonal (Figure 2a-c), square (Figure 2d-f), and isotropic array (Figure 2g-i) are compared at two different heights, h ) 10 (Figure 2a,d,g) and h ) 50 (Figure 2c,f,i). Figure 2b,e,h shows the center of the crack patterns obtained at h ) 10 zoomed by a factor of 5. This factor implies that, in the
Figure 2. Configurations at the end of the simulation of isotropic stress crack patterns (parameters: sx ) 0.1, sy ) 0.1, κ ) 0.25): the panels in the first (a-c), second (d-f), and third (g-i) rows show the results using hexagonal, square, and isotropic arrangements, respectively. Patterns for h ) 10 (panels a, b, d, e, g, and h; panels b, e, and h are 5× zooms of the centers of panels a, d, and g) and h ) 50 (panels c, f, and i) are compared.
case of a linear scaling, the average size of the crack fragments in the second and third columns should be the same. Please note that the patterns at large heights (Figure 2c,f) were already shown in our articles49,50 and are reprinted here for the sake of comparison. We have verified that the use of a larger block array does not change the average properties of the crack patterns, such as the average surface area limited by the cracks.
Simulations of Cracks Supported by Experiments
Figure 3. Configurations at the end of the simulation of unilateral stress crack patterns (parameters: sx ) 0.0, sy ) 0.1, κ ) 0.25): the panels in the first (a-c), second (d-f), and third (g-i) rows show the results using hexagonal, square, and isotropic arrangements, respectively. Patterns for h ) 10 (panels a, b, d, e, g, and h; panels b, e, and h are 5× zooms of the centers of panels a, d, and g) and h ) 50 (panels c, f, and i) are compared.
The comparison of the structures observed for hexagonal and square arrays indicates that the directions of cracks are markedly influenced by the regular arrangements of the blocks. With the isotropic grid, the crack patterns are closer to what is usually obtained from experiments with a large number of 90° angles. The comparison of the structures obtained at h ) 10 and h ) 50 shows an increase in the average surface area limited by the cracks. Moreover, the magnification of the crack patterns obtained at h ) 10 shows that the crack propagation becomes less regular and straight at smaller heights. Figure 3 shows the crack patterns observed at the end of the simulation, when the array is extended only in one direction (sx ) 0 and sy ) 0.1). The crack patterns at h ) 10 (Figure 3a,d,g) are compared to those at h ) 50 (Figure 3c,f,i) for hexagonal, square, and isotropic arrangements. Figure 3b,e,h again shows the crack patterns obtained at h ) 10 zoomed by a factor of 5. An array of, more or less linear, well-spaced cracks in the x direction appears in the hexagonal case. Please note that the stress is applied parallel to the lines of blocks building the hexagonal array. If the stress is applied perpendicular to this direction, the cracks are no longer perpendicular to the stress but form an angle of 60° with the x axis. In the case of a square array, the cracks are often not continuous, in particular, at low heights. This is explained by an inefficient propagation of the stress built up at the tip of cracks, which is due to the fact that the square array is less connected than the hexagonal one. It is remarkable that, also with the isotropic grid, linear cracks are obtained at a greater height. This shows that this kind of arrangement reproduces the external stress and not the internal arrangement of points within the grid. The comparison of the patterns obtained at h ) 10 and h ) 50 shows an increase in the average distance between the cracks. As in the case of an isotropic stress, the crack propagation becomes less regular at smaller heights. 3.2. Scaling of Crack Patterns with Height. To study the evolution of the crack patterns with height, the average area A of the fragments limited by the cracks is determined in the case of an isotropic stress. For the linear cracks produced by a unilateral stress, the patterns are analyzed using the average distance D between the cracks. From the simulations, a log-log plot of D and A as a function of κh is plotted for various κ values (κ ) 0.25, 0.5, 1,
J. Phys. Chem. C, Vol. 114, No. 41, 2010 17327
Figure 4. Logarithmic scaling of the average distance D (unilateral stress) and square root of the area A (isotropic stress) of cracks with height (parameters: κ ) 0.25, 0.5, 1.0, and 3.0). The results for square (squares), hexagonal (triangles), and isotropic (crosses) arrays are compared. The inset shows the experimental dependence of the average crack distance (D), the average square root of the areas of fragments of primary (Alp, Aip) and secondary (Als, Ais) cracks on the layer height. The average crack distance, D, on the layer height, h, designated by open squares, fits a linear curve on the logarithmic scale with log D ) 0.72 + 1.00 log h (solid lines). Solid triangles and circles correspond to the average square root of the area of primary directional and isotropic 2D cracks, Alp and Aip. The average square root of the area of the secondary directional and isotropic 2D crack domains, Als and Ais, on the layer height, h, designated by open triangles and circles, fits a linear curve on the logarithmic scale with log As ) 0.42 + 1.00 log h (dashed lines).
and 3) in Figure 4. To study the evolution over 3 orders of magnitude, a recently developed multiscale approach is applied for κh values larger than 25, as explained in the Appendix and ref 50. All data are well fitted by straight lines with a slope of 1. For the hexagonal and square arrays, their linear scaling was already observed in refs 17, 49, and 50. Here, it is shown that, in simulations using an isotropic block arrangement, there is also a linear increase in D and A with h and κ. For the square and isotropic arrays, D and A are identical within the simulation errors at a given height. Surprisingly, for the hexagonal array, A is significantly smaller than D. The large number of triangles in these crack patterns (see Figure 2a-c and the discussion in section 3.3) explains this difference because triangles have a smaller surface area than the foursided surfaces. At this stage, an important question arises: Do these simulations reflect a universal scaling law, and what are the consequences? In other words, can the scaling law observed in the simulations be proven for the model used here? Leung and Ne´da17 correctly argue that the surface area limited by cracks increases with h2, as observed in their and our simulations. Their arguments are based on how the blocks collectively slip during the simulation. We will now derive a proof of the scaling laws with height based on a new argument. To do this, five steps have to be made, as outlined in Figure 5: (i) Let us first consider a crack simulation carried out for a given height h > 1 (in Figure 5, h ) 3). The initial array of blocks is shown in Figure 5 (left grid in the first row). The simulation can be performed using reduced units, where all lengths are divided by a new scale (see step 1 in Figure 5). h is chosen as the scale. The variables in reduced units denoted by asterisks are given by the equations ri* ) ri/h, aij* ) aij/h, l* ) l/h, k* ) k, hij* ) hij, κ* ) κ. Please note that the “height” hij does not change using reduced units because it is defined as a dimensionless variable. The simulations using reduced units give
17328
J. Phys. Chem. C, Vol. 114, No. 41, 2010
the same crack patterns as those carried out without reduced units. The only difference is that all lengths are divided by the scale h. In particular, the surfaces limited by the cracks and the average crack distances are related by the equations Ah* ) Ah/ h2 and Dh* ) Dh/h. Ah* and Dh* are the properties obtained using reduced units, whereas Ah and Dh are those obtained from simulations without reduced units. (ii) The simulations can be carried out with any value for k without changing the final crack pattern. This is explained by the fact that the crack patterns obtained by the simulations depend only on the ratio of the forces Fi and Fij determined by κ and not on their absolute values. Here, we choose k*′ ) k/h (step 2 in Figure 5). The variables of the new simulation denoted by asterisks and a prime are given by the equations ri*′ ) ri/h, aij*′ ) aij/h, l*′ ) l/h, k*′ ) k/h, hij*′ ) hij, κ*′ ) κ. As discussed above, the surface areas and distances of cracks do not change on changing k, thus Ah*′ ) Ah/h2 and Dh*′ ) Dh/h. (iii) Let us now consider a simulation carried out for h ) 1 (left grid, second row in Figure 5). The simulation can be performed describing the film of particles by a finer grid of blocks (see step 3 in Figure 5). This is achieved by reducing the block-block distance and the spring length by a factor R. The variables for the simulations using a finer grid of blocks are denoted by a prime: a′ ) a/R, l′ ) l/R. To obtain the same average crack properties within the accuracy of the method, similar conditions for slipping of blocks and breaking of springs must be applied during the simulations with the original and the finer array of blocks. Therefore, the thresholds of the slipping and cracking forces and, thus, κ must be changed for the finer grid. The threshold of the slipping force is proportional to the average surface area per block, which decreases by a factor R2 (Fs′ ) Fs/R2). According to eq 1, the threshold of the spring force is proportional to the difference (a - l), which is divided by R. Therefore, the threshold of the cracking force is also reduced by R (Fc′ ) Fc/R). Combining the results for the slipping and cracking forces, we find that κ () Fc /Fs) must be increased by the factor R to obtain the same average crack properties with the finer grid. Here, we choose R ) h and obtain the following relationships for the variables of the finer grid: ri′ ) ri/h, a′ ) a/h, l′ ) l/h, k′ ) k, hij′ ) 1, κ′ ) κh. The use of a finer grid will not change the average results within the accuracy of the method. In particular, Ah)1′ ) Ah)1 and Dh)1′ ) Dh)1. Ah)1′ and Dh)1′ are the properties obtained using the finer grid, whereas Ah)1 and Dh)1 are those obtained from simulations with the original grid. (iv) The simulation of crack patterns can also be performed replacing the single spring by a bundle of β finer springs (see step 4 in Figure 5). Consequently, the spring constant must be divided by the factor β to apply the same forces in the system. In this case, the condition of block slipping is not modified. However, the spring force Fij applied to one spring is divided by the factor β. To apply the same condition for cracking (Fij > Fc), the cracking threshold Fc and, thus, κ must be reduced by the factor β. The same average crack properties within the accuracy of the method are then obtained for the simulation with one and β springs. Here, we choose β ) h and the bundle of springs is used with the finer array of blocks obtained by step 3. This leads to the following variables for the simulations, which are denoted by two primes: ri′′ ) ri/h, a′′ ) a/h, l′′ ) l/h, k′′ ) k/h, hij′′ ) hij, κ′′ ) κ′/h ) κ. In particular, Ah)1′′ ) Ah)1 and Dh)1′′ ) Dh)1. Please note that the bundle of h springs used here does not mean that we study a film with a height of h, but a film of height h ) 1 where the block-block forces are
Richardi et al.
Figure 5. Sketch to explain the proof of the scaling law in section 3.2.
described with a higher precision. Therefore, the factor κ must be reduced by a factor h. (v) A comparison of the variables used for the simulation in steps 2 and 4 shows that both sets of variables are identical, and therefore, these simulations give the same results, in particular, Ah)1′′ ) Ah*′ and Dh)1′′ ) Dh*′. Because Ah)1′′ ) Ah)1, Ah*′ ) Ah/h2, and Ah)1′′ ) Ah*′, we conclude that Ah ) h2Ah)1. In the same way, we find Dh ) hDh)1. These are the scaling laws observed in the simulations and experiments. This shows that the scaling is due to the fact that simulations at any height h > 1 are identical to one using h ) 1. The only difference is a finer array and a larger number of finer springs used at greater heights. The proof of the scaling laws implies three interesting conclusions: (i) The same arguments can be applied to any morphology of cracks. This means that the scaling law is not restricted to isotropic and linear cracks, as studied here. The same scaling should also be observed for spirals, rings, and crack patterns obtained by anisotropic stress.4,13 The proof also does not depend on the arrangement of blocks, which explains that the same linear scaling is observed for hexagonal, square, and isotropic arrays. The only condition is that κ does not change. This is fulfilled, when the force thresholds do not depend on the height. For the slipping force Fs, this is expected because the interface between the layer and its substrate should not change with height. However, it is less evident that the cracking threshold Fc does not vary for thin films. Thus, recent continuum approaches predict an evolution of Fc with height to explain the lack of cracks observed experimentally for very thin films.45,52 (ii) The proof implies that the crack pattern morphology is independent of the height. Therefore, the simulation method in its present form does not explain changes in crack patterns with height, which were observed in some experiments. For example, ref 3 finds a larger number of 120° angles and wavelike forms of cracks when the film thickness is decreased. (iii) Some changes in the simulated crack patterns with height are observed in Figures 2 and 3 and are discussed in section 3.1. With increasing height h, the simulation of crack formation is carried out with a finer grid and a larger number of springs. Therefore, only the crack patterns obtained at a greater height are reliable, whereas the structures at lower heights are modified by the low resolution of the grid. This also explains why we do
Simulations of Cracks Supported by Experiments not analyze the crack patterns obtained at low heights in the following section, because different statistics would be due only to the artifacts of the grid. The influence of grid effects on the crack patterns will be studied in detail in a future publication.56 It is, in particular, of interest for cracks in thin crystals made of monodisperse particles (such as photonic crystals57,58), where the cracks are actually influenced by the arrangements of particles. Please note that a similar proof can be given for the fact that the crack distance linearly increases with κ. 3.3. Analysis of the Morphology of Isotropic Crack Patterns. Only the theoretical crack patterns obtained at large heights (h ) 50) are taken into account because grid artifacts deteriorate the simulations at lower heights. As we explained in the last section, within our theoretical approach, the morphology of crack patterns does not change with height and, therefore, the study at one height is sufficient. To reach higher statistics at greater heights, where the number of crack fragments is small, 10 independent simulations starting from different initial random displacements of the blocks are carried out. For the hexagonal, square, and isotropic arrangements, the crack fragments have, on average, 3.97 ( 0.01, 4.0, and 3.98 ( 0.02 sides and 5.87 ( 0.01, 6.01 ( 0.04, and 5.91 ( 0.04 neighbors, respectively. Obviously, the average number of edges and neighbors of crack fragments obtained by simulations depends very slightly on the block arrangement. The values of four edges and six neighbors were explained by a recent theory proposed by Bonn et al.37 Figure 6 shows the frequencies of fragments limited by cracks as a function of their number of edges and neighbors deduced from simulations with various block arrangements for a fixed h value (h ) 50). For the square array, the number of edges (Figure 6a) is fixed at four, whereas a distribution of edge numbers is obtained for hexagonal and isotropic arrays. Furthermore, for the square array, the distribution of the number of neighbors (Figure 6b) is very narrow compared with the hexagonal and isotropic arrangements. 3.4. Comparison with Experiments. To test the predictions made by the simulations, the theoretical results are compared to experiments. We did not study the mechanical properties of our nanocrystal films. It might be interesting that another group has measured the mechanical properties of a similar CdSe nanocrystal film and showed that it can be described as a linear elastic material.59 Crack Patterns from Experiments under Isotropic and Unilateral Stresses. The fabrication and analysis of crack patterns made of maghemite nanocrystals have been explained in refs 49 and 50 and are summarized in the Supporting Information. Typical crack patterns under isotropic stress obtained by the experiments have been shown in Figure 1a,c. For the hexagonal (Figure 2a-c) and square arrays (Figure 2d-f), the comparison is limited because the directions of cracks in these simulations are too much influenced by the arrangements of the springs, while the experimental patterns were obtained in isotropic heterogeneous films due to the disordered assembly of polydisperse nanoparticles. With the isotropic grid (Figure 2g-i), the crack patterns are closer to what is obtained from experiments with a large number of 90° angles. The linear cracks obtained by the experiments under unilateral stress (Figure 1b,d) are reproduced by the simulations for all three grids at sufficiently large heights (Figure 2c,f,i). The cracks perpendicular to the primary linear cracks observed in thick films (Figure 1d) were recently reproduced in simulations by the application of an anisotropic stress.49 Please note that such
J. Phys. Chem. C, Vol. 114, No. 41, 2010 17329
Figure 6. Average frequencies of crack fragments as a function of the number of sides (panel a) and neighbors (panel b). The results for isotropic crack patterns (parameters: sx ) 0.1, sy ) 0.1, κ ) 0.25, and h ) 50) are shown. The simulation results obtained using hexagonal, square, and isotropic arrangements of blocks are compared to the experimental values.
ladder-like crack patterns widely appear in paintings (e.g., the Mona Lisa of Leonardo da Vinci) and in nature. Experimental Scaling of Crack Distances with Height. In refs 49 and 50 (see also the Supporting Information), a linear scaling of the crack distance with height was observed, in excellent agreement with the simulations (Figure 4). This linear scaling is explained by the proof proposed in section 3.2. Obviously, simulations using the square or isotropic arrays better describe the experimental data, which show an identity of D and A.50 Therefore, these simulations are used for the quantitative comparison with experiments. The simulations for square and isotropic arrays yield an average value of 1.6 for the ratio A/(Hκ) ) D/(Hκ). Using the experimental value of K ) A/H ) D/H ) 5.4 obtained for primary cracks,50 we find a value of 3.4 for the ratio κ of the cracking to slipping force thresholds, which applies to the experiments. This κ value is compared to that deduced from a simple van der Waals model: The threshold for the cracking force is approximated by the negative adhesion force of two flat van der Waals surfaces60
Fc ) -Fad )
A 4γ ) 3 D 6πD
(2)
where A and D are the Hamaker constant and the separation between the two surfaces. The interfacial tension γ() Α/24πD2)
17330
J. Phys. Chem. C, Vol. 114, No. 41, 2010
Richardi et al.
corresponds to the change in the free energy due to the formation of the two free surfaces of the cracks. The threshold of the slipping force is calculated from the cobble stone model.61 assuming that, due to the slipping motion, the separation between the film and the substrate increases by a small amount ∆D, while the move in lateral distance is ∆σ. This leads to the threshold of the slipping, which corresponds to the friction force61
Fs )
4γ′ε∆D D′∆σ
(3)
where the interfacial tension γ′ is the free energy change due to the separation of the nanocrystal film and the substrate necessary for the slipping. ε is the fraction of the above surface energy, which is lost every time the film moves across the characteristic length ∆σ. Combining eqs 2 and 3 yields
κ)
Fc γD′∆σ ) Fs Dγ′ε∆D
(4)
Assuming typical values for the parameters γ ≈ γ′ (≈ 25 mJ m-2), D ≈ D′ (≈ 0.2 nm), ∆D ≈ 0.05 nm, ∆σ ≈ 0.1 nm, and ε ) 0.1-0.5,60 the κ values are found in the range from 4 to 20. The κ value of 3.4 obtained by the comparison of experiment and theory is at the lower limit of this range, which may be explained by an underestimation of the interfacial tension γ′, which describes the contact between the massive substrate and the hydrocarbons coating the nanocrystals. Indeed, in our model, it is estimated that this tension γ′ is the same as γ, which describes the contact between two nonmassive van der Waals walls made of the hydrocarbons coating the assembled nanocrystals. However, it can be expected that the interfacial tension for a massive substrate (γ′) is larger than that for a particle aggregate (γ), which leads to a decrease in the estimated value for κ, according to eq 4, in good agreement with our results. Analysis of the Morphology of Experimental Isotropic Crack Patterns. Finally, let us consider the details of the morphology of isotropic crack patterns (e.g., the numbers of sides and neighbors) produced from the experiments. Figure 7a-c shows the angle distributions of the crack junctions at three different film heights: 2.8, 5.9, and 18.4 µm. The angle distribution is always peaked at 90°. In very thin films (Figure 7a), the distribution becomes asymmetrical with a larger probability of finding junctions with angles up to 120°. From a height of 5.9 µm, the angle distribution is independent of the film thickness with the central peak (90° ( 2.5°) between 70 and 80%. Figure 7d-f shows the distributions of the numbers of sides (edges) and neighbors at the same heights. For a height of 2.8 µm (Figure 7d), note that a rather large fraction of five-sided fragments exists (35%) in coincidence with a relatively high fraction of five neighbors (24%). On increasing the film height from 2.8 µm (Figure 7e) to 18.4 µm (Figure 7f), the percentage of four-sided fragments markedly increases (from 53% to 80%), whereas the number of five-sided fragments drops (35% to 13%) and the number of three-sided fragments remains quite constant (between 12% and 7%). The relative percentage of six neighbors decreases from 42% to 38%, the numbers of fragments with five and seven neighbors remain quite constant (around 25% and 18%, respectively). Above a height of 18.4 µm, the distributions of edges and neighbors no longer change within the statistical accuracy of the results. In Figure 6, the averaged frequencies for heights above 5.9 µm are plotted. Note that the
Figure 7. Distribution of crack junction angles and n-sided polygons in the final isotropic crack patterns for selected samples. The average layer heights are 2.8 ( 0.02 (a, d), 5.9 ( 0.06 (b, e), and 18.4 ( 0.10 µm (c, f). Filled and open bars show the frequencies of sides and neighbors.
experimental values are in good agreement with the results of Bohn et al.37 measured for cracks in ceramic glazes. Shorlin et al.4 observed a slightly lower frequency of about 60% for foursided fragments in his experiments. It is interesting to note that this value coincides with our values observed in very thin films. The experimental results are now compared with the simulations. The experiments have shown that the average numbers of edges and neighbors of crack fragments do not depend on the film thickness. Average experimental values of 4.08 ( 0.03 sides and 5.94 ( 0.05 neighbors are obtained, in good agreement with the simulation results. Figure 6 shows that the square array predicts frequencies for both the edges and the neighbors that markedly deviate from the experimental results. From Figure 6b, a rather good agreement between simulated data obtained from hexagonal and isotropic grids with experiments is observed. However, the side number distribution with hexagonal grids (Figure 6a) is rather broad compared with the experimental results, indicating that the best array to describe the experimental crack pattern is the isotropic grid. This shows that the isotropy within layers markedly influences the frequency of edges and neighbors of crack fragments. With respect to the dependence of the crack statistics on the experimental parameters, it usually does not depend on the height except for very thin films. It should also not depend on the solvent and the coating of the particle because both were changed without any effect on the statistics. With respect to the theoretical approach, it is important to use a sufficient number of springs to avoid grid effects. For the same reason, very small values of κ must be avoided. Please note also that the polydispersity used to obtain the isotropic grid must be larger than 7% to avoid ordered self-assembly of the particles constituting the grid. 4. Conclusions We have used simulations to systematically study the formation of crack patterns in drying films. The thickness of
Simulations of Cracks Supported by Experiments the film is varied by a factor up to 1000, and the anisotropy of the patterns is changed by applying unilateral and isotropic stresses. Although the simulations reproduce the increase of crack distances with height observed in experiments, the crack propagation is markedly influenced by the spring arrangements used in the simulations. This has no impact on the average number of neighbors and sides of the crack fragments obtained by simulations using an isotropic stress. Thus, the fragments have, on average, four edges and six neighbors, in good agreement with the experiment, which is expected for hierarchical structures with 90° junctions. However, the frequency of fragments as a function of the number of edges or neighbors is markedly influenced by the arrangements of blocks. Only simulations using an isotropic grid yield correct results in comparison with experiments: This shows that a correct treatment of the isotropy is essential to obtaining a qualitative and quantitative reproduction of crack patterns. Experimentally, the same linear scaling law with height is found in the isotropic and directional cases. This is correctly reproduced using square or isotropic arrangements, but not for the hexagonal one. A simple van der Waals model can explain the ratio of the slipping and cracking threshold, κ, deduced from a comparison of the experiments and the simulations. We have derived a proof for the linear scaling law based on the use of reduced units. This shows that, within our approach, there is no difference between simulations carried out at different heights except for a different resolution of the cracks due to a finer array of blocks. This proof implies that the scaling of crack patterns is really a universal feature valid for any morphology of the cracks. Within our model, a change in morphology of the crack patterns cannot be induced by a variation of the height. The difference in the cracks observed between the simulations at various heights is due to the different resolution of the grid with respect to the size of the cracks. Actually, in the experiments, the average geometries of crack patterns observed in the experiments, such as the angles of junctions and the number of sides of crack fragments, do not depend on the height except within very thin layers. Appendix Details of the Simulation Method. Generation of Isotropic Grids. The generation of isotropic arrangements placing the blocks in a random manner is usually accompanied by the appearance of holes within the grid that would be nucleation centers for cracks. To avoid these holes, we carry out a Monte Carlo simulation of 1000 particles having various sizes in two dimensions using periodic boundary conditions. The particles interact by a very weak Lennard-Jones (LJ) potential (interaction parameters: σ j ) 1.0, ε ) 0.1 in reduced units). A Gaussian distribution with a variation of 10% describes the size distribution of the particle diameter. A size distribution is used to avoid an ordered hexagonal assembly of the particles, which would give an anisotropic hexagonal grid. Starting from an initial random position of particles, the length of the box is reduced every 10 000 steps by a factor 0.99. At the end of the simulation run, a highly compressed system of particles is obtained with an average distance of 1.0. The centers of particles are duplicated in both directions to obtain an array of box length 400. Each pair of neighboring centers i and j is connected by h springs whenever their distance is smaller than 1.5 (σi + σj), where σi and σj are the diameters of particles i and j. Thus, an isotropic array of points without predominant orientations is obtained. Because of its high compacity, the holes within the grids are avoided. We investigated the size distribution between 0 and
J. Phys. Chem. C, Vol. 114, No. 41, 2010 17331 30% and observe that the simulation results do not depend on the size distribution when it is sufficiently large (>8%), as presented here (10%). Details of the Simulations. A simulation step consists of the following: (1) It starts with the calculation of the forces of all blocks in the array. (2) Both thresholds, Fs and Fc, are lowered, keeping κ constant until either threshold is exceeded somewhere in the system. (3) If Fi > Fs the block i is slipped to a new position where Fi ) 0 using a Newton-Raphson method.62 (4) If Fij > Fc, a spring is broken (hij f hij - 1). (5) The forces of the neighboring blocks are corrected and blocks are moved or springs are broken until all forces in the system respect the conditions Fi < Fs and Fij < Fc. To avoid the instantaneous breaking of the whole bundle of springs before the neighboring stress has at all relaxed, only one spring of a bundle is broken during a simulation step. Each simulation takes 10 million steps (several hours on a modern workstation). In the study of isotropic crack patterns, about 10 000 blocks (Lx ) Ly ) 100) are sufficient to determine the scaling law. However, for the simulations of linear 1D cracks, the number of blocks must be drastically increased. This is explained as follows: The number of surfaces limited by cracks for isotropic patterns is proportional to the square of the side length of the array. In contrast, the number of bands limited by linear cracks in the case of unilateral stress increases only linearly with the side length. For example, for κ ) 0.25 and h ) 40 using 10 000 blocks, an isotropic crack pattern contains about 11 surfaces, whereas the linear crack pattern consists only of about three bands. Therefore, to reach similar statistical accuracies in the 1D and 2D cases, the size of the array has to be increased. Here, an array of about 250 000 blocks (Lx ) Ly ) 500; the exact number of blocks depends on the extension) is used in comparison to 10 000 blocks for 2D crack patterns used in ref 17. It should be noted that, because of a larger array, the number of simulation steps must also be increased from 300 000 to 10 million steps. Because of the increased size of the array and the large step number, the simulation time greatly increases. The most time is spent by the program in searching the largest values of Fi and Fij in the array (step 2). Producing lists containing the blocks with the largest values of Fi and Fij can considerably enhance this step. This list must be updated only for the forces that change their values in steps 3-5. Thus, the simulation time is reduced by 2 orders of magnitude (still several hours on a modern workstation). Leung and Ne´da17 use an approximation to speed up the force calculations in the case of an isotropic stress (see eq 1 in ref 17). To test this approximation, simulations are carried out with and without the approximation. It turns out that, for an isotropic extension, the approximation does not influence the simulation results. However, for a unilateral or anisotropic stress, the approximated forces are very different from the exact ones, leading to deviations in the final crack patterns. This is due to the fact that the approximation is based on an isotropic expansion using the positions of the initial array, which is not correct in the case of an anisotropic stress. All simulations were carried out using the homemade C package, CrackT. Superblock Method. In ref 41, a new method, the so-called superblock approach, was developed to achieve simulations varying the height by several orders of magnitude. Within the superblock approach, the distance and spring constant between
17332
J. Phys. Chem. C, Vol. 114, No. 41, 2010
super blocks in comparison to the original blocks are increased by a factor denoted by R > 1. To obtain the same average crack patterns at given values of h and κ for the superblocks as for the original block array, κ must be multiplied by the factor R, as explained in ref 49. Acknowledgment. The authors gratefully acknowledge the contribution of Henning Kayser to this work. We also like to thank Emmanuel Trizac and Jacob N. Israelachvili for fruitful discussions. Supporting Information Available: The experimental method, including chemicals and apparatus, synthesis of solutions of maghemite nanocrystals, fabrication and analysis of crack patterns, and the scaling of the average crack distance in experiments. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Walker, J. Sci. Am. 1986, 225, 178. (2) Herrmann, H. J.; Roux, S. Statistical Models for the Fracture of Disordered Media; Elsevier Science: Amsterdam, 1990. (3) Groisman, A.; Kaplan, E. Europhys. Lett. 1994, 25, 415. (4) Shorlin, K. A.; de Bruyn, J. R.; Graham, M.; Morris, S. W. Phys. ReV. E 2000, 61, 6950. (5) Pauchard, L.; Lazarus, V.; Abou, B.; Sekimoto, K.; Aitken, G.; Lahanier, C. Reflets Phys. 2007, 3, 5. (6) Ruiz-Agudo, E.; Martin-Ramos, J. D.; Rodriguez-Navarro, C. J. Phys. Chem. B 2006, 110, 3727–3733. (7) Teoh, M. M.; Chung, T. S.; Pramoda, K. P. J. Phys. Chem. B 2006, 110, 5889. (8) Lisiecki, I.; Albouy, P. A.; Pileni, M. P. J. Phys. Chem. B 2004, 108, 20050. (9) Zhao, X. S.; Siepmann, J. I.; Xu, W.; Kiang, Y. H.; Sheth, A. R.; Karaborni, S. J. Phys. Chem. B 2009, 113, 5929. (10) Georges, S.; Ayres Rocha, R.; Djurado, E. J. Phys. Chem. C 2008, 112, 3194. (11) Hu, Z.; Chen, Q.; Li, Z.; Yu, Y.; Peng, L. M. J. Phys. Chem. C 2010, 114, 881. (12) Lintanf, A.; Mantoux, A.; Blanquet, E.; Djurado, E. J. Phys. Chem. C 2007, 111, 5708. (13) Allain, C.; Limat, L. Phys. ReV. Lett. 1995, 74, 2981. (14) Leung, K. T.; Jozsa, L.; Ravasz, M.; Ne´da, Z. Nature 2001, 410, 166. (15) Li, H. L.; Dong, W.; Bongard, H. J.; Marlow, F. J. Phys. Chem. B 2005, 109, 9939–9945. (16) Komatsu, T. S.; Sasa, S. Jpn. J. Appl. Phys. 1997, 36, 391. (17) Leung, K. T.; Ne´da, Z. Phys. ReV. Lett. 2000, 85, 662. (18) Meakin, P. Thin Solid Films 1987, 151, 165. (19) Skjeltorp, A. T.; Meakin, P. Nature 1982, 335, 424. (20) Andersen, J. V.; Brechet, Y.; Jensen, H. J. Europhys. Lett. 1994, 26, 13. (21) Andersen, J. V. Phys. ReV. B 1994, 49, 9981. (22) Leung, K. T.; Andersen, J. V. Europhys. Lett. 1997, 38, 589. (23) Hornig, T.; Sokolov, I. M.; Blumen, A. Phys. ReV. E 1996, 54, 4293. (24) Crosby, K. M.; Bradley, R. M. Phys. ReV. E 1997, 55, 6084. (25) Walmann, T.; Malthe-Sørenssen, A.; Feder, J.; Jøsang, T.; Meakin, P.; Hardy, H. H. Phys. ReV. Lett. 1996, 77, 5393.
Richardi et al. (26) Malthe-Sørenssen, A.; Walmann, T.; Feder, J.; Jøssang, T.; Meakin, P.; Hardy, H. H. Phys. ReV. E 1998, 58, 5548. (27) Ne´da, Z.; Leung, K. T.; Jo´zsa, L.; Ravasz, M. Phys. ReV. Lett. 2002, 88, 095502. (28) Hayakawa, Y. Phys. ReV. E 1994, 49, R1804. (29) Hayakawa, Y. Phys. ReV. E 1994, 50, R1748. (30) Marder, M. Phys. ReV. E 1994, 49, R51. (31) Sasa, S.; Sekimoto, K.; Nakanishi, H. Phys. ReV. E 1994, 50, R1733. (32) Yuse, A.; Sano, M. Nature 1993, 362, 329. (33) Colina, H.; de Arcangelis, L.; Roux, S. Phys. ReV. B 1993, 48, 3666. (34) Morgenstern, O.; Sokolow, I. M.; Blumen, A. Europhys. Lett. 1993, 22, 487. (35) Morgenstern, O.; Sokolow, I. M.; Blumen, A. J. Phys. A 1993, 26, 4521. (36) Sokolov, I. M.; Morgenstern, O.; Blumen, A. Macromol. Symp. 1994, 81, 235. (37) Bohn, S.; Douady, S.; Courder, Y. Phys. ReV. Lett. 2005, 94, 054503. (38) Bohn, S.; Platkiewicz, J.; Etreotti, B.; Adda-Bedia, M.; Couder, Y. Phys. ReV. E 2005, 71, 046215. (39) Bohn, S.; Pauchard, L.; Couder, Y. Phys. ReV. E 2005, 71, 046214. (40) Sherman, D.; Be’ery, I. Phys. ReV. Lett. 2004, 93, 265501. (41) Livne, A.; Ben-David, O.; Fineberg, J. Phys. ReV. Lett. 2007, 98, 124301. (42) Kierfeld, J.; Vinokur, V. M. Phys. ReV. Lett. 2006, 96, 175502. (43) Pilipenko, D.; Spatschek, R.; Brener, E. A.; Mu¨ller-Krumbhaar, H. Phys. ReV. Lett. 2007, 98, 015503. (44) Måløy, K. J.; Hansen, A.; Hinrichsen, E. L.; Roux, S. Phys. ReV. Lett. 1992, 68, 213. (45) Singh, K. B.; Tirumkudulu, M. S. Phys. ReV. Lett. 2007, 98, 218302. (46) Bouchaud, E.; Lapasset, G.; Plane`s, J. Europhys. Lett. 1990, 13, 73. (47) Termonia, Y.; Meakin, P. Nature 1986, 320, 429. (48) Vashishta, P.; Kalia, R. K.; Nakano, A. J. Phys. Chem. B 2006, 110, 3727–3733. (49) Ngo, A. T.; Richardi, J.; Pileni, M. P. Nano Lett. 2008, 8, 2485. (50) Ngo, A. T.; Richardi, J.; Pileni, M. P. J. Phys. Chem. B 2008, 112, 14409. (51) Pileni, M. P.; Lalatonne, Y.; Ingert, D.; Lisiecki, I.; Courty, A. Faraday Discuss. 2004, 125, 251. (52) Jia, S.; Banerjee, S.; Herman, I. P. J. Phys. Chem. C 2008, 112, 162–171. (53) Wang, X.; Qi, Z.; Zhao, C.; Wang, W.; Zhang, Z. J. Phys. Chem. C 2009, 113, 13139–13150. (54) Zhang, X.; Imae, T. J. Phys. Chem. C 2009, 113, 5947–5951. (55) Pauchard, L.; Elias, F.; Boltenhagen, P.; Cebers, A.; Bacri, J. C. Phys. ReV. E 2008, 77, 021402. (56) Richardi, J.; Ngo, A. T.; Pileni, M. P. Manuscript in preparation. (57) Freeman, D.; Madden, S.; Luther-Davies, B. Opt. Express 2005, 13, 3079. (58) Scharrer, M.; Wu, X.; Yamilov, A.; Cao, H.; Changa, R. P. H. Appl. Phys. Lett. 2005, 86, 151113. (59) Banerjee, S.; Jia, S.; Kim, S. I.; Robinson, R. D.; Kysar, J. W.; Bevk, J.; Herman, I. P. Nano Lett. 2006, 6, 175. (60) Hunter, R. J. Foundations of Colloid Science; Oxford University Press: Oxford, U.K., 1986. (61) Israelachvili, J. Surface Forces and Nanorheology of Molecular Thin Films. In Springer Handbook of Nanotechnology, 2nd ed.; Bhushan, B., Ed.; Springer: Heidelberg, Germany, 2007; pp 1-66. (62) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; Cambridge University Press: New York, 2000.
JP104329G
