Drying Patterns: Sensitivity to Residual Stresses

0.71 MB·English

by  Yossi Cohen

Drying Patterns: Sensitivity to Residual Stresses
Yossi Cohen1, Joachim Mathiesen2, and Itamar Procaccia1
1Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel
2Physics of Geological Processes, University of Oslo, Oslo, Norway
(Dated: January 7, 2009)
Volume alteration in solid materials is a common cause of material failure. Here we investigate
thecrack formation in thin elastic layers attached toasubstrate. Weshowthat small variations in
thevolumecontraction and substraterestraint can producewidely differentcrack patternsranging
fromspiralstocomplexhierarchicalnetworks. Thenetworksareformedwhenthereisnoprevailing
9
gradient in material contraction whereas spirals are formed in the presence of a radial gradient in
0
thecontraction of a thin elastic layer.
0
2

n Introduction. Desiccationis known to produce com- the substrate. In this process the coating often shrinks
a plex networks of shrinkage-cracks in starch-water mix- and tensile stresses are produced [7] and if one is less
J
tures or clays[1, 2, 3, 4]. In concrete small cracks are careful, the stress may result in an unwanted cracking
7
often formed by the preparatory drying process and by of the coating. Due to the spinning of the system, the

the later ingress of reactive reagents. Similarly in na- residual stress of the coating may also contain inherent
]
i ture, the infiltration of fluids and chemical reagents into shear components. Other mechanisms like anisotropic
c
s rocks generate internal stresses that form intricate pat- drying rates can also leave behind remnant shear. In
-
l terns of pervasive cracks[5]. Typically the stress is gen- caseswherethethicknessofthedryingspecimenissmall
r
eratedfromlocalvolumechanges. Fracturesare alsoob- and the contractionfairly uniform, i.e. no residualshear
t
m
servedin thin films attached to a substrate. Experiment stresses,thegrowingcrackstypicallyformanintricatehi-
. on films have revealed intricate patterns ranging from erarchicalpattern. The patternis the resultofa cascade
t
a the hierarchicalstructure typically observed in mud and of successive cracks (which is supported by experimen-
m
concrete to spiral shaped cracks [1, 6]. In spin-coating talevidences[2]); ateachfragmentationstage,a crackis
- a fluid droplet is added at the center of a rotating sub- forming which divides a mother fragment of area A into
d
strate and is spread by centrifugal forces to cover the two daughter fragments of areasA and A respectively,
n 1 2
o full substrate. During the drying and curing of the sys- with areaconservation(A=A1+A2). It is worthnotic-
c tem,chemicalbondsareformedbetweenthecoatingand ing that in principle the trajectory of the crack that is
[ dividingthemotherfragmentcanbeanything,andisde-

1 termined from the shape and size of the mother domain
v and the inherent material disorder.
7
In order to analyze the spiral and hierarchical cracks,
9
we consider a system consisting of a thin elastic layer
7
0 attachedto anelastic substrate. Under plane stress con-
. ditions we have that the in-plane strain tensor ǫ in the
1 ij
0 layer is related to the stress tensor σij by
9
0
ǫ = (σ νσ )/E+β ,
: xx xx− yy
v
ǫ = (σ νσ )/E+β ,
yy yy xx
i −
X ǫ = (1+ν)σ /E, (1)
xy xy
r
a
where E is Young’s modulus and β (in the absence of
external stress) is a measure of the free volume change
caused by e.g. drying or thermal expansion of the thin
elastic layer. Whenever the film is displaced from its
equilibrium 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-
FIG.1: Coloronline. Fourdifferentstagesintheevolutionof
early proportional to u. In general it is assumed that
a hierarchical crack network. The total contraction was 9%
−
andcrackswerenucleatedinsidedomainswheneverthemax- volume alteration in the film happens on a time scale
imum principal stress exceeded σc = 0.85 in arbitrary units. much larger than the time required for elastic waves to
The substrate restraining force was drawn from a uniform propagate across the system and the system is therefore
distribution with 10% disorder and a mean of unity. assumed to always be in elastostatic equilibrium. The
2
force balance therefore assumes the form nowprovideanestimateofacriticaldomainsizethatwill
fracture under a predefined yield stress. As long as this
∂jσij µui =0 , (2) yieldstress is lowerthanthe materialstressfracture will
−
form and grow. Depending on the material contraction
whereµistheconstantofproportionalityofthesubstrate
the fractures may develop into spiral shaped patterns or
restoringforce. Forsmalldeformationsthestrainfollows
hierarchicalnetworks. Inbothcasesthemaximumstress
from the displacement via the relations ǫ = (∂ u +
ij j i is reduced by the propagating crack and only when the
∂ u )/2. Combining this relation with the force balance
i j stressdropsbelowtheyieldstressthefracturingstops. In
Eq. (2)andstress-strainrelationsEq. (1)weachievethe
Fig. 1weshowafracturenetworkresultingfromnumer-
following equation for the displacement
ical solutions as explained below, formed from an initial
contractionofthe substrateandapredefinedyieldstress
1+ν 2(1+ν)µ
u+ ( u)= u . (3) level. AccordingtoEq. (7)eachdomaindivisionreduces
△ 1 ν∇ ∇· E
− the stress and an average linear size R of the domains
h i
We now provide an estimate of the typical stress en- can be found for a given yield stress by inverting Eq.(7).
countered during volume alteration of thin films with a The non-uniformity of the elastic layer and the com-
linearspatialextendofsizeR. Tothatend,weshallcon- plex boundary conditions make it hard andoften impos-
siderthe maximumstressfora circulardomainofradius sible to find an analytical solution to the displacement
Rlocatedatthecenterofcoordinatesandwithvanishing equation. Therefore we have implemented a numerical
stress at the boundaries. The displacement field is for a method based on the Galerkin finite element discretiza-
uniformmaterialcontractionβ foundasasolutiontothe tionusinganadaptivetriangularmeshing. Inthevicinity
radial symmetric version of Eq. (3) of a propagating crack tip we highly increase the resolu-
tion by decreasing locally the area of the triangular ele-
∂2ur 1∂ur 1 ments and thereby allow for an accurate computation of
+ a+ u =0 , (4)
∂r2 r ∂r −(cid:18) r2(cid:19) r thestressintensityfactorsofthepropagatingcrack. The
drying process is simulated by applying a body force to
where the material specific constant a is given by a = the elements, i.e. we shift the equilibrium position by
(1 ν2)µ/E. MultiplyingbothsidesofEq. (4)byr2 and adding an extra force term on the right hand side of Eq.
−
rescalingr with√ayieldsthemodifiedbesseldifferential (2). Inthatwaywecanreadilyadddisorderintothesys-
equation. The solution for σrr(R) = 0 and ur(0) = 0 is tem by selecting the magnitude of the local body force
given by from a random distribution. In the simulations on hi-
erarchical fracture networks presented below, we use a
βRI (√ar)
u (r)= 1 , u =0. uniform distribution with unit mean.
r θ
√aRI0(√aR) (1 ν)I1(√aR) Crack Initiation and propagation. Here we
− −
(5)
presentindetailhow wenucleate cracksandmodeltheir
Here I , n = 0,1 is the modified bessel function of the
n evolution. First we find the points which have the high-
first kind. From the displacement field, the stress in
est stress and exceed the critical value. The stress is
cylindrical coordinates follows from the expressions
determined along the principal axes of the stress matrix
σ , i.e. the principal stress. Whenever the yield stress
E ∂u u ij
σrr(r) = 1 ν2 (cid:18)∂r +νr −β(1+ν)(cid:19), is exceeded, we nucleate at the point of yielding a small
− semi elliptical void with an eccentricity of 0.998. The
E u ∂u
majoraxis of the ellipse is alignedin the directionof the
σ (r) = +ν β(1+ν) . (6)
θθ 1 ν2 (cid:18)r ∂r − (cid:19) maximum principal stress. We allow the crack to evolve
−
according to the Griffith criterion and the principle of
Note that by the symmetry of the problem the shear
local symmetry, i.e. the crack will grow in a direction
stress vanishes. The breaking of this symmetry will be
such as to annul the local shear component at the crack
important for the formation of the spiral crack patterns
tip. At each step of propagation we compute the stress
presentedbelow. Thestresscomponentshavetheirmax-
in every element near the crack tip and find the stress
imum (absolute value) at the middle of the circular do-
intensity factors from a best fit to the equations [14].
main and are given by
K θ K θ θ
σ = I cos3 3 II sin cos2 , (8)
σ (0)=σ (0)= βE 1+ν 1 . (7) θθ √2πr 2 − √2πr 2 2
rr θθ
1−ν (cid:18)2I0(√aR) − (cid:19) σ = KI sinθ cos2 θ + KII cosθ(1 3sin2 θ) .
rθ
√2πr 2 2 √2πr 2 − 2
The magnitude of the stress components monotonically
increaseswithRandinthelimitR the stresscom- Here r,θ are local polar coordinates with respect to the
→∞
ponents achieve the value Eβ/(1 ν). In the limit crack tip with θ measured from the line following the
− −
R=0 the stress becomes Eβ/2. From Eq. (7) we can directionofthecrack. σ andσ arethecircumferential
θθ rθ
−
3
Contraction
0
0 0.
20.0 )distribution −1.0 9% 17% 23%
2 | ) 005.00 (ogCum. 10−2.0 100
1 2. L 0
3.
−
c-P( | 0.50 −2.4 L−o2g.100(Area)−1.6 ation 120 0
t
0 o
2 R
0. 5%
10%
05 15% 3
0.
100
0.0 0.1 0.2 0.3 0.4
| c- 1 2 |
FIG. 2: Color online. The distribution density of |χ−1/2|
for simulations with 5%,10% and 15% disorder on the sub- FIG. 3: Simulation of spiral cracks for various values of the
strate restraining. The distributions are averaged over do- material contraction β and the rotation angles. In all the
mains formed at the 6th generation of cracks. No significant panels, thecrack was initiated at the centerand was allowed
variation is seen at these fairly low levels of disorder. The to propagate until it reached the outer boundary. Note that
black line on top represents a best fit with an exponential the smaller the shear stress is the more pronounced is the
distribution exp(−|χ−1/2|/α) where α=0.03. In the inset spiralling.
weshowforthesamedatathecumulativedistributionsofthe
domain areas. The dashed line on top is an estimate of the
distributions considering the individual domain divisions to the maximum of the strain energy release rate
be uncorrelated (see text). dG(α)/dα=0 is equivalent to K (α)=0 or dK (α)/dα=
II I
0, thus the new direction of the crack α corresponds
0
to the point where K (α ) exhibits a maximum and
tensile stress and the shear stress, respectively. K and I 0
I
K (α )=0[13]. Applying the latterto Eq. (10)yields,
K are the unknown stress intensity factors for mode II 0
II
I and II, respectively. The principle of local symmetry
is satisfied if the crack grows in a direction given by an α0 =2arctan(cid:18)(KI −qKI2+8KI2I)/4KII(cid:19) . (13)
angle α where K 0. Suppose that the crack forms
II
→
aninfinitesimalkinkatanangleαfromtheolddirection We emphasize that in this model both the cracking time
of the crack,we can define the localmode I and mode II and the area of the fragmented elements depend only
stress intensity factors, on material contraction and initial disorder. In a natu-
ral system this may not always be the case since mate-
K (α) = limσ √2πr (9)
I r→0 θθ rial properties and disorder can evolve in time. Uniform
θ θ θ contractionof the elastic layer produces homogenous hi-
= K cos3 3K sin cos2 .
I II erarchical crack patterns. One can in this case argue
2 − 2 2
K (α) = limσ √2πr (10) that the effect of the crack is to partition the mother
II rθ
r→0 area A (of generic shape) into two areas A1 = χA and
θ θ θ θ A = (1 χ)A, where 0 < χ < 1 is a random variable
= KIsin2cos2 2 +KIIcos2(1−3sin2 2) . wh2ose dis−tribution (that must be symmetric under the
Whether the crack propagates or not is dictated by the transformation χ 1 χ) is unknown. In Fig. 2 is the
7→ −
Griffith criterion, i.e. the energy balance of the energy distribution of χ 1/2 = A1 A2 /A shown together
| − | | − |
releaserateintothecracktipregionmustbalancethedis- with a best fit to an exponential distribution. Although
sipation involvedin the crack propagation. For a kinked thedomainareasarecorrelatedtotheirmotherdomains,
crack, the energy release rate is, the exponential distribution of χ allows for a simple es-
timate of the areadistribution by neglecting the correla-
G(α)=(K2(α)+K2 (α))/E . (11)
I II tion. Thatistheareasatthen’thgenerationlevelarecan
Since bedeterminedbyaproductofnrandomnumbersdrawn
dKI= 3K cos2 αsinα 3K cosα(1 3sin2 α)= 3K ,fromthe exponentialdistribution,i.e. Ai(n) =A0 nj χij.
dα −2 I 2 2−2 II 2 − 2 −2 II In Fig. 2, we show in the inset a distribution forQdomain
(12) areasatthe6thgenerationtogetherwiththedistribution
4
5 elastic layer relative to the underlying substrate with an
0.
angle θ (r) = θ rγ), the crack would propagate along
9% eq 0
a spiral trajectory. Different powers of γ result in the
13%
0.4 17% formation of different spirals. For the simulation of the
20%
crack propagationwe use an initially circular symmetric
system. Acrackis then initiated atthe centerof the cir-
3
0. cle and is allowed to propagate according to the Griffith
criterionuntilitreachestheboundary. Theresultsusing
k
arotationofthesubstratebyapowerγ =1/2areshown
2
0. 0.06 inFig. 3usingvariousvaluesfortheprefactorsθ0 andβ,
respectively. The cracks have a shape that fit well a log-
3
1 kb0.0 arithmic spiral, i.e. they have a form r(θ) = Aexp(κθ)
0.
0 where κ depends on the material contraction β and the
0
0.
0 1 2 3 4 5 6 7 rotation θ . In Fig. 4 we show best fits of κ as function
Rotation angle [p 100] 0
0.0 of θ0 and for four values of β.
0 1 2 3 4 5 6 7
Insummary,wepointedouttheroleofresidualstresses
Rotation angle [p 100]
in determining the crack patterns in drying thin sub-
strates. For spiral patterns we related the properties of
FIG.4: Coloronline. Thefigureshowstherelationbetweenκ
the spiral to the degree of residual shear stress left in
andtherotationangleforfourdifferentvaluesofthematerial
the layer. For hierarchical patterns we determined the
contraction β. κ is computedas theexponentof a best fit to
a logarithmic spiral, r(φ)=r0eκφ. The rotation angle is the position where new cracks initiate as a function of the
prefactor θ0 used in the expression for the relative rotation mother cell, and offered a relation of final mean size of
between the substrate and the thin film. The inset shows a cells to the critical value of the yield stress.
data collapse of κβ for the same curves and is in agreement
with thesimple scaling form κ∼θ0/β. We thank M. Adda-Bedia for proposing the study of
spiralcracks. ThisworkhasbeensupportedbythePGP,
a Center of Excellence at the University of Oslo, the
oftheproductofrandomnumbers. Afterafewnumberof
GermanIsraeliFoundationandtheMinervaFoundation,
generations this distribution will approach a log-normal
Munich, Germany.
distribution. The deviation from the exponential distri-
butionforlargervalues of χ 1/2 willforanincreasing
| − |
number of fracture generations lead to a less good fit
using the exponential distribution as an approximation.
Spirals. We now investigate what happens when the [1] A. T. Skjeltorp and P. Meakin, Nature 335, 424 (1988).
contraction is non-uniform and has smooth gradients. [2] S. Bohn, J. Platkiewicz, B. Andreotti, M. Adda-Bedia,
and Y.Couder, Phys.Rev.R71, 046215 (2005)
Experiments on thin films attached to a substrate by a
[3] S. Bohn, S. Douady and Y. Couder, Phys. Rev. Lett.,
spin-coatingtechniquerevealabroadrangeofcrackpat-
94, 054503 (2005).
terns[6]rangingfromnetworksofcracksto singlecracks [4] A. Groisman and E. Kaplan, Europhys. Lett., 25, 415
spiraling outwardsfromtheir site of nucleation. The nu-
(1994).
cleationisusuallytakingplaceatlocalizedsiteswithhigh [5] K. Iyer, B. Jamtveit, J. Mathiesen, A. Malthe-Srenssen
stresstypicallygeneratedfromsmalldefectsorinclusions and J. Feder, EPSL 267, 503 (2008).
in the material. If the material contraction is uniform [6] M.SendovaandK.Willis,Appl.phys.A,76,957(2003).
in the neighborhood around the crack, the crack would [7] J. Malzbender and G. de With, Thin Solid Films 359,
210 (2000).
propagatestraighttowardsthe materialboundarywhere
[8] E. Katzav, M. Adda-Bedia and B. Derrida, Europhys.
it would curve to meet the boundary at a right angle. Lett., 78, 46006 (2007).
However if the material contraction increases smoothly [9] S.Sadhukhan,S.R.Majumder,D.Mal,T.DuttaandS.
away from the site of nucleation, straight cracks would Tarafdar, J. Phys. Cond. Matt., 19, 356206 (2007).
be unstable and wouldstartto curve. In that waycircu- [10] J.V.AndersenandL.J.Lewis,Phys.Rev.E57,R1211
larshapedcrackscanbeformed. Duringthepreparation (1998)
of thin film coatings (such as spin-coating) it is not un- [11] A. Buchel and J. P. Sethna, Phys. Rev. E 55, R7669
(1997)
commontohaveaminorresidualshearstress. Theshear
[12] J. Liang, R. Huang, J. H. Prevost and Z. Suo, Int. J.
stressbreaksthesymmetryofthesystemandthecircular Solids Struct.40, 2343 (2003).
crack may then turn into a spiral. If we alter the con- [13] D. Broek, Elementary engineering fracture mechanics.
tractionsuchthatitincreaseslinearlyawayfromagiven Kluwer Academic Publishers, Dordrecht,1986.
siteofcracknucleation,(e.g. followsasimplelinearform [14] K. B. Broberg, Crack and Fracture, (Academic Press,
β(r) = βr, and add a small shear stress by rotating the London) 1999.