Drying Patterns: Sensitivity to Residual Stresses
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n Drying Patterns: Sensitivity to Residual Stresses
Yossi Cohen , Joachim Mathiesen , and Itamar Procaccia Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel Physics of Geological Processes, University of Oslo, Oslo, Norway (Dated: November 23, 2018)Volume alteration in solid materials is a common cause of material failure. Here we investigatethe crack formation in thin elastic layers attached to a substrate. We show that small variations inthe volume contraction and substrate restraint can produce widely different crack patterns rangingfrom spirals to complex hierarchical networks. The networks are formed when there is no prevailinggradient in material contraction whereas spirals are formed in the presence of a radial gradient inthe contraction of a thin elastic layer.
Introduction . Desiccation is known to produce com-plex networks of shrinkage-cracks in starch-water mix-tures or clays[1, 2, 3, 4]. In concrete small cracks areoften formed by the preparatory drying process and bythe later ingress of reactive reagents. Similarly in na-ture, the infiltration of fluids and chemical reagents intorocks generate internal stresses that form intricate pat-terns of pervasive cracks[5]. Typically the stress is gen-erated from local volume changes. Fractures are also ob-served in thin films attached to a substrate. Experimenton films have revealed intricate patterns ranging fromthe hierarchical structure typically observed in mud andconcrete to spiral shaped cracks [1, 6]. In spin-coatinga fluid droplet is added at the center of a rotating sub-strate and is spread by centrifugal forces to cover thefull substrate. During the drying and curing of the sys-tem, chemical bonds are formed between the coating and
FIG. 1: Color online. Four different stages in the evolution ofa hierarchical crack network. The total contraction was 9%and cracks were nucleated inside domains whenever the max-imum principal stress exceeded σ c = 0 .
85 in arbitrary units.The substrate restraining force was drawn from a uniformdistribution with 10% disorder and a mean of unity. the substrate. In this process the coating often shrinksand tensile stresses are produced [7] and if one is lesscareful, the stress may result in an unwanted crackingof the coating. Due to the spinning of the system, theresidual stress of the coating may also contain inherentshear components. Other mechanisms like anisotropicdrying rates can also leave behind remnant shear. Incases where the thickness of the drying specimen is smalland the contraction fairly uniform, i.e. no residual shearstresses, the growing cracks typically form an intricate hi-erarchical pattern. The pattern is the result of a cascadeof successive cracks (which is supported by experimen-tal evidences [2]); at each fragmentation stage, a crack isforming which divides a mother fragment of area A intotwo daughter fragments of areas A and A respectively,with area conservation ( A = A + A ). It is worth notic-ing that in principle the trajectory of the crack that isdividing the mother fragment can be anything, and is de-termined from the shape and size of the mother domainand the inherent material disorder.In order to analyze the spiral and hierarchical cracks,we consider a system consisting of a thin elastic layerattached to an elastic substrate. Under plane stress con-ditions we have that the in-plane strain tensor ǫ ij in thelayer is related to the stress tensor σ ij by ǫ xx = ( σ xx − νσ yy ) /E + β ,ǫ yy = ( σ yy − νσ xx ) /E + β ,ǫ xy = (1 + ν ) σ xy /E, (1)where E is Young’s modulus and β (in the absence ofexternal stress) is a measure of the free volume changecaused by e.g. drying or thermal expansion of the thinelastic layer. Whenever the film is displaced from itsequilibrium position by a local displacement u the elas-tic substrate tries to restore the film by a force f ( u ).For small displacements we assume that this force is lin-early proportional to − u . In general it is assumed thatvolume alteration in the film happens on a time scalemuch larger than the time required for elastic waves topropagate across the system and the system is thereforeassumed to always be in elastostatic equilibrium. Theforce balance therefore assumes the form ∂ j σ ij − µu i = 0 , (2)where µ is the constant of proportionality of the substraterestoring force. For small deformations the strain followsfrom the displacement via the relations ǫ ij = ( ∂ j u i + ∂ i u j ) /
2. Combining this relation with the force balanceEq. (2) and stress-strain relations Eq. (1) we achieve thefollowing equation for the displacement △ u + 1 + ν − ν ∇ ( ∇ · u ) = 2(1 + ν ) µE u . (3)We now provide an estimate of the typical stress en-countered during volume alteration of thin films with alinear spatial extend of size R . To that end, we shall con-sider the maximum stress for a circular domain of radius R located at the center of coordinates and with vanishingstress at the boundaries. The displacement field is for auniform material contraction β found as a solution to theradial symmetric version of Eq. (3) ∂ u r ∂r + 1 r ∂u r ∂r − (cid:18) a + 1 r (cid:19) u r = 0 , (4)where the material specific constant a is given by a =(1 − ν ) µ/E . Multiplying both sides of Eq. (4) by r andrescaling r with √ a yields the modified bessel differentialequation. The solution for σ rr ( R ) = 0 and u r (0) = 0 isgiven by u r ( r ) = βRI ( √ ar ) √ aRI ( √ aR ) − (1 − ν ) I ( √ aR ) , u θ = 0 . (5)Here I n , n = 0 , σ rr ( r ) = E − ν (cid:18) ∂u∂r + ν ur − β (1 + ν ) (cid:19) ,σ θθ ( r ) = E − ν (cid:18) ur + ν ∂u∂r − β (1 + ν ) (cid:19) . (6)Note that by the symmetry of the problem the shearstress vanishes. The breaking of this symmetry will beimportant for the formation of the spiral crack patternspresented below. The stress components have their max-imum (absolute value) at the middle of the circular do-main and are given by σ rr (0) = σ θθ (0) = βE − ν (cid:18) ν I ( √ aR ) − (cid:19) . (7)The magnitude of the stress components monotonicallyincreases with R and in the limit R → ∞ the stress com-ponents achieve the value − Eβ/ (1 − ν ). In the limit R = 0 the stress becomes − Eβ/
2. From Eq. (7) we can now provide an estimate of a critical domain size that willfracture under a predefined yield stress. As long as thisyield stress is lower than the material stress fracture willform and grow. Depending on the material contractionthe fractures may develop into spiral shaped patterns orhierarchical networks. In both cases the maximum stressis reduced by the propagating crack and only when thestress drops below the yield stress the fracturing stops. InFig. 1 we show a fracture network resulting from numer-ical solutions as explained below, formed from an initialcontraction of the substrate and a predefined yield stresslevel. According to Eq. (7) each domain division reducesthe stress and an average linear size h R i of the domainscan be found for a given yield stress by inverting Eq.(7).The non-uniformity of the elastic layer and the com-plex boundary conditions make it hard and often impos-sible to find an analytical solution to the displacementequation. Therefore we have implemented a numericalmethod based on the Galerkin finite element discretiza-tion using an adaptive triangular meshing. In the vicinityof a propagating crack tip we highly increase the resolu-tion by decreasing locally the area of the triangular ele-ments and thereby allow for an accurate computation ofthe stress intensity factors of the propagating crack. Thedrying process is simulated by applying a body force tothe elements, i.e. we shift the equilibrium position byadding an extra force term on the right hand side of Eq.(2). In that way we can readily add disorder into the sys-tem by selecting the magnitude of the local body forcefrom a random distribution. In the simulations on hi-erarchical fracture networks presented below, we use auniform distribution with unit mean. Crack Initiation and propagation . Here wepresent in detail how we nucleate cracks and model theirevolution. First we find the points which have the high-est stress and exceed the critical value. The stress isdetermined along the principal axes of the stress matrix σ ij , i.e. the principal stress. Whenever the yield stressis exceeded, we nucleate at the point of yielding a smallsemi elliptical void with an eccentricity of 0 . stressintensity factors from a best fit to the equations [14]. σ θθ = K I √ πr cos θ − K II √ πr sin θ θ , (8) σ rθ = K I √ πr sin θ θ K II √ πr cos θ − θ . Here r, θ are local polar coordinates with respect to thecrack tip with θ measured from the line following thedirection of the crack. σ θθ and σ rθ are the circumferential . . . . . . | c - P ( | c- | ) −2.4 −2.0 −1.6 − . − . − . . Log ( Area ) Log ( C u m . d i s t r i bu t i on ) FIG. 2: Color online. The distribution density of | χ − / | for simulations with 5% ,
10% and 15% disorder on the sub-strate restraining. The distributions are averaged over do-mains formed at the 6th generation of cracks. No significantvariation is seen at these fairly low levels of disorder. Theblack line on top represents a best fit with an exponentialdistribution exp( −| χ − / | /α ) where α = 0 .
03. In the insetwe show for the same data the cumulative distributions of thedomain areas. The dashed line on top is an estimate of thedistributions considering the individual domain divisions tobe uncorrelated (see text). tensile stress and the shear stress, respectively. K I and K II are the unknown stress intensity factors for modeI and II, respectively. The principle of local symmetryis satisfied if the crack grows in a direction given by anangle α where K II →
0. Suppose that the crack formsan infinitesimal kink at an angle α from the old directionof the crack, we can define the local mode I and mode IIstress intensity factors, K I ( α ) = lim r → σ θθ √ πr (9)= K I cos θ − K II sin θ θ .K II ( α ) = lim r → σ rθ √ πr (10)= K I sin θ θ K II cos θ − θ . Whether the crack propagates or not is dictated by theGriffith criterion, i.e. the energy balance of the energyrelease rate into the crack tip region must balance the dis-sipation involved in the crack propagation. For a kinkedcrack, the energy release rate is, G ( α ) = ( K I ( α ) + K II ( α )) /E . (11)Since dK I dα = − K I cos α α − K II cos α − α − K II , (12) %9 %17 %23 Contraction R o t a t i o n FIG. 3: Simulation of spiral cracks for various values of thematerial contraction β and the rotation angles. In all thepanels, the crack was initiated at the center and was allowedto propagate until it reached the outer boundary. Note thatthe smaller the shear stress is the more pronounced is thespiralling. the maximum of the strain energy release rate dG ( α ) /dα =0 is equivalent to K II ( α )=0 or dK I ( α ) /dα =0, thus the new direction of the crack α correspondsto the point where K I ( α ) exhibits a maximum and K II ( α ) = 0 [13]. Applying the latter to Eq. (10) yields, α = 2 arctan (cid:18) ( K I − q K I + 8 K II ) / K II (cid:19) . (13)We emphasize that in this model both the cracking timeand the area of the fragmented elements depend onlyon material contraction and initial disorder. In a natu-ral system this may not always be the case since mate-rial properties and disorder can evolve in time. Uniformcontraction of the elastic layer produces homogenous hi-erarchical crack patterns. One can in this case arguethat the effect of the crack is to partition the motherarea A (of generic shape) into two areas A = χA and A = (1 − χ ) A , where 0 < χ < χ − χ ) is unknown. In Fig. 2 is thedistribution of | χ − / | = | A − A | /A shown togetherwith a best fit to an exponential distribution. Althoughthe domain areas are correlated to their mother domains,the exponential distribution of χ allows for a simple es-timate of the area distribution by neglecting the correla-tion. That is the areas at the n’th generation level are canbe determined by a product of n random numbers drawnfrom the exponential distribution, i.e. A ( n ) i = A Q nj χ ij .In Fig. 2, we show in the inset a distribution for domainareas at the 6th generation together with the distribution . . . . . . Rotation angle [ p k . . . Rotation angle [ p kb FIG. 4: Color online. The figure shows the relation between κ and the rotation angle for four different values of the materialcontraction β . κ is computed as the exponent of a best fit toa logarithmic spiral, r ( φ ) = r e κφ . The rotation angle is theprefactor θ used in the expression for the relative rotationbetween the substrate and the thin film. The inset shows adata collapse of κβ for the same curves and is in agreementwith the simple scaling form κ ∼ θ /β . of the product of random numbers. After a few number ofgenerations this distribution will approach a log-normaldistribution. The deviation from the exponential distri-bution for larger values of | χ − / | will for an increasingnumber of fracture generations lead to a less good fitusing the exponential distribution as an approximation. Spirals . We now investigate what happens when thecontraction is non-uniform and has smooth gradients.Experiments on thin films attached to a substrate by aspin-coating technique reveal a broad range of crack pat-terns [6] ranging from networks of cracks to single cracksspiraling outwards from their site of nucleation. The nu-cleation is usually taking place at localized sites with highstress typically generated from small defects or inclusionsin the material. If the material contraction is uniformin the neighborhood around the crack, the crack wouldpropagate straight towards the material boundary whereit would curve to meet the boundary at a right angle.However if the material contraction increases smoothlyaway from the site of nucleation, straight cracks wouldbe unstable and would start to curve. In that way circu-lar shaped cracks can be formed. During the preparationof thin film coatings (such as spin-coating) it is not un-common to have a minor residual shear stress. The shearstress breaks the symmetry of the system and the circularcrack may then turn into a spiral. If we alter the con-traction such that it increases linearly away from a givensite of crack nucleation, (e.g. follows a simple linear form β ( r ) = βr , and add a small shear stress by rotating the elastic layer relative to the underlying substrate with anangle θ eq ( r ) = θ r γ ), the crack would propagate alonga spiral trajectory. Different powers of γ result in theformation of different spirals. For the simulation of thecrack propagation we use an initially circular symmetricsystem. A crack is then initiated at the center of the cir-cle and is allowed to propagate according to the Griffithcriterion until it reaches the boundary. The results usinga rotation of the substrate by a power γ = 1 / θ and β ,respectively. The cracks have a shape that fit well a log-arithmic spiral, i.e. they have a form r ( θ ) = A exp( κθ )where κ depends on the material contraction β and therotation θ . In Fig. 4 we show best fits of κ as functionof θ and for four values of β .In summary, we pointed out the role of residual stressesin determining the crack patterns in drying thin sub-strates. For spiral patterns we related the properties ofthe spiral to the degree of residual shear stress left inthe layer. For hierarchical patterns we determined theposition where new cracks initiate as a function of themother cell, and offered a relation of final mean size ofcells to the critical value of the yield stress.We thank M. Adda-Bedia for proposing the study ofspiral cracks. This work has been supported by the PGP ,a Center of Excellence at the University of Oslo, theGerman Israeli Foundation and the Minerva Foundation,Munich, Germany. [1] A. T. Skjeltorp and P. Meakin, Nature , 424 (1988).[2] S. Bohn, J. Platkiewicz, B. Andreotti, M. Adda-Bedia,and Y. Couder, Phys.Rev.R , 046215 (2005)[3] S. Bohn, S. Douady and Y. Couder, Phys. Rev. Lett., , 054503 (2005).[4] A. Groisman and E. Kaplan, Europhys. Lett., , 415(1994).[5] K. Iyer, B. Jamtveit, J. Mathiesen, A. Malthe-Srenssenand J. Feder, EPSL , 503 (2008).[6] M. Sendova and K. Willis, Appl. phys. A, , 957 (2003).[7] J. Malzbender and G. de With, Thin Solid Films ,210 (2000).[8] E. Katzav, M. Adda-Bedia and B. Derrida, Europhys.Lett., , 46006 (2007).[9] S. Sadhukhan, S. R. Majumder, D. Mal, T. Dutta and S.Tarafdar, J. Phys. Cond. Matt., , 356206 (2007).[10] J. V. Andersen and L. J. Lewis, Phys. Rev. E , R1211(1998)[11] A. Buchel and J. P. Sethna, Phys. Rev. E , R7669(1997)[12] J. Liang, R. Huang, J. H. Prevost and Z. Suo, Int. J.Solids Struct. , 2343 (2003).[13] D. Broek, Elementary engineering fracture mechanics .Kluwer Academic Publishers, Dordrecht, 1986.[14] K. B. Broberg,