WIXLIB

ticity, modeled by a mechanical analog consisting of a spring anda dashpot in parallel, for which the relaxation modulus is 共t兲 共t兲共18兲where is the stiffness of the spring, representing the elasticshear modulus at the rubbery limit, and is the viscosity. TheLaplace transform of the relaxation modulus is 共s兲 s共19兲After substituting 共19兲 into Eqs. 共15兲 and 共16兲, inverse Laplacetransform leads to1 4 H2 S3 u1 H u1 S1 t 4共1 兲 x1 1 4 H2 S1 1 2 H u3S3 u3 t 2共1 兲 4共1 兲 x1 coupled equations, 共20兲 and 共22兲, must be used in this case. Generalization of the plane-strain response to the three-dimensionalwould take the similar form as the lubrication theory 关14兴, but willnot be further pursued in the present study.2.3 Coupled Evolution Equations. The interface betweenthe elastic and viscoelastic layers is assumed to remain bondedduring deformation. Consequently, the displacements and tractions are continuous across the interface, which couples the equilibrium equations of the elastic layer, Eqs. 共1兲 and 共2兲, with thetime-dependent responses of the viscoelastic layer, Eqs. 共24兲 and共25兲, and leads to冉1 2 H 4w 2w w N w Df N x x x x x x x x t 2共1 兲 共20兲 共21兲 w 冊共26兲Similarly, for an incompressible viscoelastic layer 共 0.5兲, theinverse transform of Eq. 共17兲 leads to u H N u t x u3H 3 2S 3 H 2 S 1 u3 t3 x21 2 x1 Equations 共26兲 and 共27兲 are coupled, nonlinear evolution equations, which may be solved numerically to simulate threedimensional deformation of an elastic-viscoelastic bilayer andevolution of the resulting two-dimensional wrinkle patterns. In theremainder of this paper, however, we focus our attention on planestrain deformation and one-dimensional wrinkles only, leaving thetwo-dimensional wrinkles for a subsequent study. The reducedequations for the plane-strain wrinkles are summarized as follows:共22兲Equations 共20兲 and 共22兲 have the similar form as the Reynold’slubrication theory for nearly parallel flow of a thin liquid layer关14,20兴, but with an additional term accounting for the elasticlimit of the viscoelastic layer.For the present study, we assume a compressible viscoelasticlayer 共i.e., 0.5兲 and further neglect the H2 terms in Eqs. 共20兲and 共21兲 for thin-layer approximation, which leads to u1 H u , S1 t 1共23兲1 2 H u3 S3 u t 2共1 兲 3共24兲Here the two traction components are assumed to have comparable magnitudes and the thickness of the viscoelastic layer isassumed to be small compared to the wavelength 共L 2 / k兲.Equation 共23兲 is equivalent to a shear lag model, which assumesuniform shear stress across the thin layer. Similar models havebeen used for both elastic and viscous layers 关21,22兴. Equation共24兲 is similar to the Winkler model for elastic foundation 关23兴 butincludes a time derivative term due to the viscous effect. The twoequations are uncoupled under the thin-layer approximation.In the above development, plane-strain deformation and periodic surface tractions have been assumed. The restriction of periodic tractions has been relaxed by using differentiation of thesurface tractions with respect to x1 in places of the particular wavenumber after inverse Laplace transform. The resulting equations,共20兲 and 共21兲, are apparently independent of wave number andcan be used for arbitrary tractions by linear superposition of theirFourier components. At the end, the in-plane and out-of-planeresponses are decoupled by the thin-layer approximation. Therefore, the restriction of plane-strain deformation can be relaxed bygeneralizing the in-plane response, Eq. 共23兲, for both x1 and x2directions, which leads to u H u S t Journal of Applied Mechanics冉冊1 2 H 4w 2w N w w Df 4 N 2 w 共28兲 x x t 2共1 兲 x x u H N u t x N 0h f Efhf1 2f共29兲冋 冉 冊册 u 1 w x 2 x2共30兲Recently, Huang et al. 关24兴 developed a similar model to simulate the evolution of two-dimensional wrinkle patterns in elasticfilms on soft substrates, where the viscoelastic Kelvin model wasused to relate the lateral deflection and the normal traction, similarto Eq. 共24兲, but the relation between the in-plane displacement andthe shear traction was taken to be linear elastic. While the attention there was focused on various equilibrium wrinkle patterns,the interest of the present study is the temporal evolution ofwrinkles.3Linear Perturbation AnalysisAssume a small deflection of the elastic layer in the form of共31兲w共x,t兲 A共t兲cos kxFor the linear perturbation analysis, the evolution of the in-planedisplacement is uncoupled from the lateral deflection and, therefore, ignored. Inserting 共31兲 into Eq. 共28兲 and retaining only theleading-order terms in A, we obtain that共25兲for 1 , 2. Equations 共24兲 and 共25兲 then represent the approximate solution for the three-dimensional response of a thin viscoelastic layer subjected to the boundary conditions in Eqs. 共6兲and 共7兲.In the case of an incompressible viscoelastic layer, however, theWinkler-type equation for the out-of-plane displacement 共Eq.共24兲兲 breaks down and the decoupling is not applicable. The共27兲dA E f A共t兲 dt where 共1 2 兲k2Hh f24共1 兲共1 2f 兲冋 k2h2f 12共1 2f 兲 0Ef共32兲册共33兲Solving Eq. 共32兲 leads toNOVEMBER 2005, Vol. 72 / 957

lengths. On the other hand, the fastest growing wavelength doesnot change, but the corresponding growth rate decreases. The fastest growth rate reduces to zero at a critical ratio冉 冊 EfFig. 2 Initial growth rate as a function of wavelength by thelinear perturbation analysis, for various ratios between the rubbery modulus of the viscoelastic layer and the Young’s modulus of the elastic layer冉冊A共t兲 A0 exp st 共34兲where A0 is the initial perturbation amplitude, / E f is a characteristic time scale, and s / E f is the dimensionless growthrate of the perturbation. The stability of the bilayer therefore depends on the sign of the growth rate. If s 0 for all wave numbersk, the bilayer is stable and remains flat. Otherwise, if s 0 for anypermissible wave numbers, the bilayer is unstable and perturbations grow to form wrinkles. In this case, the amplitude growsexponentially with time at the initial stage. A more sophisticatedanalysis 关16兴 has shown that the initial growth can be nonexponential if the viscoelastic layer has a finite elastic modulus at theglassy state 共elastic limit as t 0兲.Figure 2 plots the growth rate as a function of the perturbationwavelength 共L 2 / k兲 for different ratios between the rubberymodulus of the viscoelastic layer and the elastic modulus of theelastic layer. At the limiting case when 0, s and thegrowth rate is positive 共recall that 0 0兲 for long wave perturbations, as shown by the dashed line in Fig. 2. Consequently, thebilayer is always unstable. The critical wavelength isLc h f冑 Ef3共1 2f 兲 0共35兲which is identical to the critical length of Euler buckling. Thegrowth rate is positive when L Lc and peaks at the wavelengthLm h f冑 2E f3共1 2f 兲 0共36兲Similar results were obtained for an elastic film on a viscous layer关14,25兴, where the fastest growing wavelength Lm is, however,shorter by 13.4% due to the incompressibility of the viscous layer.Using typical values for a thin aluminum layer: E f 70 GPa, f 0.35, h f 40 nm, and 0 100 MPa, we obtain Lc 2.05 mand Lm 2.90 m. The latter compares closely to the initial wavelengths observed in experiments by Yoo and Lee 关15兴 despite therough estimate of the stress.As the ratio / E f increases, the growth rate decreases; thecurve in Fig. 2 shifts down, but without any change in the shape.As a result, the critical wavelength increases, and a second criticalwavelength emerges at the long wave end. The growth rate is nowpositive within a window bounded by the two critical wave958 / Vol. 72, NOVEMBER 2005 c冉 冊3共1 2f 兲共1 2 兲 H 02共1 兲hf Ef2共37兲The bilayer becomes stable when / E f is greater than the criticalratio. Alternatively, Eq. 共37兲 may be rewritten to give the criticalcompressive stress, below which a bilayer with the given thickness ratio and moduli ratio is stable. The critical condition is identical to that for an elastic film on a thin elastic substrate with theshear modulus 关12,16兴.It is noted that, by the critical condition in Eq. 共37兲, the stabilityof an elastic-viscoelastic bilayer depends on the rubbery modulus共i.e., the long-term limit of the relaxation modulus兲 of the viscoelastic layer, but independent of the initial modulus 共e.g., theglassy state兲. In other words, despite that the viscoelastic layer isinitially stiff or even rigid, the bilayer “foresees” the subsequentsoftening of the layer and becomes unstable spontaneously. Thetime scale of wrinkle growth is proportional to the viscosity, andthe growth rate increases as the rubbery modulus decreases. Thewavelength of the fastest growing mode, however, is independentof the viscoelastic layer, as given in Eq. 共36兲. Our previous study关16兴 showed that the fastest growing wavelength weakly dependson the thickness ratio and Poisson’s ratio. The thin-layer approximation in the present study leads to a reasonably accurate wavelength, but underestimates the growth rate for the fastest growingmode when the thickness ratio H / h f is larger than 2.4Equilibrium WrinklesSetting / t 0 in Eqs. 共28兲 and 共29兲 leads to two coupled nonlinear ordinary differential equations, from which one can solvefor equilibrium states. The solution is identical to that for an elastic film on a thin elastic substrate with the shear modulus .The latter has been obtained by an energy minimization procedure关12,16兴, as summarized below. First, the equilibrium amplitude ofa sinusoidal wrinkle with a wave number k is given byAeq 冋2冑1 2f2共1 兲 1 0共kh f 兲2 kE f 12共1 2f 兲1 2 E f k2Hh f册1/2共38兲It can be confirmed that only when the bilayer is unstable doesthere exist nonzero, real-valued equilibrium wrinkle amplitudes.Furthermore, minimization of the elastic strain energy in the bilayer with respect to the wave number selects an equilibriumwrinkle wavelengthLeq h f冋2共1 2 兲Ef H2 3共1 兲共1 f 兲 h f册1/4共39兲The corresponding wrinkle amplitude can be obtained from Eq.共38兲 with k 2 / Leq. Again, using typical values: E f 70 GPa, f 0.35, h f 40 nm, 0.01 MPa, 0.45, H 400 nm, and 0 100 MPa, we obtain Leq 7.00 m and Aeq 71.9 nm. It isnoted that Eq. 共39兲 underestimates the equilibrium wavelengthwhen the thickness ratio H / h f is larger than 2.Comparing the equilibrium wrinkle wavelength to the initiallyfastest growing wavelength given by Eq. 共36兲, we note that thetwo wavelengths can be totally independent. The fastest growingwavelength, which dominates the initial growth, is determined bythe kinetics and depends on the compressive stress in the elasticlayer, but independent of the viscoelastic layer. The equilibriumwavelength, on the other hand, is determined by energetics anddepends on the thickness and rubbery modulus of the viscoelasticlayer, but independent of the stress in the elastic layer. Such independence may enable simultaneous determination of the residualstress and rubbery modulus from the initial and final wrinklewavelengths, respectively.Transactions of the ASME

Fig. 3 Evolution of the lateral deflection w and the in-planedisplacement u by numerical simulation with a sinusoidal initial perturbationAt the equilibrium state, the shear traction at the interface isnearly zero and the lateral displacement approximately takes theform 关12,14兴1u 8 kA2eq sin共2kx兲共40兲where k 2 / Leq. The wavelength of the in-plane displacement ishalf of the wrinkle wavelength at the equilibrium state.Fig. 4 Amplitude of a sinusoidal wrinkle as a function of timeJournal of Applied MechanicsFig. 5 Numerical simulation of evolving wrinkles with a random initial perturbation. The left column shows the deflectionof the elastic layer, and the right column shows the corresponding Fourier spectra.5Numerical SimulationsIn this section, we simulate the evolution of wrinkles by numerically integrating the nonlinear equations, 共28兲 and 共29兲. Forsimplicity, we use the explicit forward-time-center-space 共FTCS兲finite difference method. The algorithm is conditionally stable. Toachieve sufficient accuracy, a small space step x is first specified.Next, the time step t is determined by the stability and convergence of the numerical results. In the following simulations, weuse x 1.0h f and t 0.1 / E f . The time is normalized by thetime scale / E f , which ranges widely from 1 s to 1 s, depending on the material of the viscoelastic layer and the temperature. In all simulations, the periodic boundary condition is assumed.The bilayer is in equilibrium at the reference state 共Fig. 1共a兲兲with no tractions at the interface. By introducing a small perturbation displacement to the reference state as the initial condition,the system evolves until it reaches another equilibrium state. First,we start with a sinusoidal deflection of amplitude A0 0.01h f andzero in-plane displacement at t 0. The wavelength L 30h f wasselected to be close to the fastest growing wavelength 共Lm 26.9h f 兲 to save the computational time. Other parameters are 0 0.01E f , 0.0001E f , H 10h f , f 0.3, and 0.45. Asnoted before, the present model underestimates the growth rateand the equilibrium wavelength for thick viscoelastic layers 共H 2h f 兲. Nevertheless, the wrinkling kinetics should be similar, andNOVEMBER 2005, Vol. 72 / 959

Fig. 6 Evolution of the dominant wrinkle wavelength by numerical simulation: „a short time evolution and „b long timeevolutionin the numerical simulations, we use H 10h f as in the metal/polymer bilayer experiments by Yoo and Lee 关15兴. Figure 3 showssnapshots of the evolving displacements. The amplitude of thelateral deflection grows with time, but the wavelength remainsconstant for the entire period of the simulation up to t 50,000 .Meanwhile, relatively small in-plane displacement evolves concomitantly, but with a wavelength half of the wrinkle wavelength,as predicted by the equilibrium solution in Eq. 共40兲. Figure 4shows the wrinkle amplitude as a function of time. The amplitudefirst grows exponentially, as predicted by the linear perturbationanalysis 共shown as the straight dashed line兲. Starting at about t 20,000 , the growth rate deviates from the linear behavior andgradually approaches a plateau. The amplitude essentially remainsconstant after t 40,000 , indicating that an equilibrium state hasbeen reached. The equilibrium amplitude given by Eq. 共38兲 for theselected wavelength 共L 30h f 兲 is Aeq 0.537h f , as indicated by thehorizontal dotted line in Fig. 4. The result from the numericalsimulation agrees closely with the analytical solutions by the linear perturbation analysis at the initial stage and by the energeticanalysis for the equilibrium amplitude.In the above simulation, the wrinkle wavelength is arbitrarilyselected a priori and the evolution stops when the correspondingequilibrium state is reached. However, the wavelength is not necessarily the equilibrium wavelength selected by energy minimization, as given in Eq. 共39兲. In other words, the equilibrium statereached in the previous simulation is energetically unstable. Tofurther relax the strain energy, continual evolution is possible oncethe equilibrium state is perturbed with different wavelengths 关14兴.In real situations, various sources 共e.g., thermal fluctuation and960 / Vol. 72, NOVEMBER 2005Fig. 7 The root mean square „RMS of the wrinkle as a function of time: „a short-time evolution and „b long-timeevolutionsurface defects兲 may induce the initial perturbation, which is random, in general. Figure 5 shows a numerical simulation of evolving wrinkles that starts from a random initial perturbation. Thein-plane displacement is again zero initially 共not shown兲. Otherparameters are 0 0.01E f , 0.00001E f , H 10h f , f 0.3,and 0.45. For each snapshot of the wrinkle, the correspondingFourier spectrum is shown to the right. Although many wavelengths coexist in the initial perturbation, only those of intermediate wavelengths grow and the fastest growing wavelength dominates at the initial stage. Consequently, an increasingly regularwrinkle emerges from the initially random perturbation. As theevolution continues, the amplitude grows and the wavelength increases. After a sufficiently long time, only one wavelength remains and the wrinkle reaches an equilibrium state. Figure 6 plotsthe evolution of the dominant wavelength 共maximum intensity inthe Fourier spectrum兲, and Fig. 7 shows the root-mean square共rms兲 of the wrinkle as a function of time. Also plotted in Figs. 6and 7 are the simulation results with a larger rubbery modulus, 0.0001E f , for comparison.From the numerical simulations, three stages of wrinkle evolution can be identified: initial growth of the fastest growing mode,intermediate growth with mode transition, and, finally, an equilibrium wrinkle state. Such behavior qualitatively agrees with theexperimental observations in a metal/polymer bilayer film 关15兴. Atthe initial stage, the wavelength of the fastest growing mode predicted by the linear perturbation analysis is Lm 26.9h f , which isindependent of the rubbery modulus. Figure 6共a兲 shows that thedominant wavelengths in the two simulations are indistinguishable up to t 2 104 , and the wavelength is close to the predictedvalue. During this stage, the wrinkle amplitude grows exponenTransactions of the ASME

tially, but the growth rate depends on the rubbery modulus. InFigure 7共a兲, the two dashed lines indicate the exponential growthpredicted by the linear perturbation analysis, where the largermodulus leads to slower growth. At the intermediate stage, boththe amplitude and wavelength of the wrinkle evolve toward theequilibrium state. The analytical solutions for the equilibrium stateare indicated as dashed lines in Figs. 6共b兲 and 7共b兲. For 0.00001E f , the equilibrium wrinkle has a wavelength Leq 60.0h f and an amplitude Aeq 1.63h f 共rms 1.15h f 兲. For 0.0001E f , the equilibrium wrinkle has a wavelength Leq 33.7h f and an amplitude Aeq 0.619h f 共rms 0.438h f 兲. The equilibrium states agree closely with the numerical results. It is notedthat, although the initial growth is slower, the time to reach theequilibrium state is significantly shorter with the larger rubberymodulus for the viscoelastic layer.6SummaryWe have developed a nonlinear model for temporal evolution ofwrinkles in elastic-viscoelastic bilayer thin films. The modelcouples a nonlinear theory of elastic plates with a thin-layer approximation of linear viscoelastic responses. Although the modelis three-dimensional in nature, the analyses and numerical simulations of the present study have focused on plane-strain deformation. Analytical solutions are obtained for the linear perturbationanalysis and the equilibrium state. Numerical simulations illustrate the evolution process from the initial growth to the equilibrium state. The results show that the kinetics of wrinkling stronglydepend on the viscoelastic layer.AcknowledgmentThe authors are grateful for the support by NSF Grant No.CMS-0412851 and by the Texas Advanced Technology Program.References关1兴 Iacopi, F., Brongersma, S. H., and Maex, K., 2003, “Compressive Stress Relaxation Through Buckling of a Low-k Polymer-Thin Cap Layer System,”Appl. Phys. Lett., 8

been assumed to be a constant independent of time, considering the factor that the Poisson's ratio is typically a weak function of time. If the viscoelastic layer is incompressible i.e., 0.5 , how-ever, Eq. 16 takes a different form u 3 x 1,s 2 1 2ks s 2 3 3kH B s kH A s cos kx 1 17 where the first term in the brackets scales .