Asymptotic regimes in elastohydrodynamic and stochastic leveling on a viscous film
Christian Pedersen John Niven Thomas Salez Kari Dalnoki-Veress Andreas Carlson


Mechanics Division, Department of Mathematics, University of Oslo, 0316 Oslo, Norway.

Department of Physics and Astronomy, McMaster University, 1280 Main Street West, Hamilton, Ontario, L8S 4M1, Canada.

Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33405 Talence, France.
Global Station for Soft Matter, Global Institution for Collaborative Research and Education, Hokkaido University, Sapporo, Hokkaido 060-0808, Japan.

Department of Physics and Astronomy, McMaster University, 1280 Main Street West, Hamilton, Ontario, L8S 4M1, Canada.

acarlson@math.uio.no
Mechanics Division, Department of Mathematics, University of Oslo, 0316 Oslo, Norway.



An elastic sheet that deforms near a solid substrate in a viscous fluid is a situation relevant to various dynamical processes in biology, geophysics and engineering. Here, we study the relaxation dynamics of an elastic plate resting on a thin viscous film that is supported by a solid substrate. By combining scaling analysis, numerical simulations and experiments, we identify asymptotic regimes for the elastohydrodynamic leveling of a surface perturbation of the form of a bump, when the flow is driven by either the elastic bending of the plate or thermal fluctuations. In both cases, two distinct regimes are identified when the bump height is either much larger or much smaller than the thickness of the pre-wetted viscous film. Our analysis reveals a distinct crossover between the similarity exponents with the ratio of the perturbation height to the film height.

Asymptotic regimes in elastohydrodynamic and stochastic leveling on a viscous film
    Andreas Carlson
    December 30, 2023
==================================================================================





§ INTRODUCTION

The motion of an elastic sheet supported by a thin layer of viscous fluid is a phenomenon that manifests itself in processes spanning wide ranges of time and length scales, from e.g. magmatic intrusion in the Earth's crust <cit.>, to fracturing and crack formation in glaciers <cit.>, to pumping in the digestive and arterial systems <cit.>, or the construction of 2D crystals for electronic engineering <cit.>. Elastohydrodynamic flows have been studied in model geometries in order to understand their generic features and the inherent coupling between the driving force from the elastic deformations of the material and the viscous friction force resisting motion <cit.>.

The investigation of an initially flat elastic membrane that is subsequently subjected to an applied deformation has helped disclose how system size, magnitude and direction of elastic deformations and spatial confinement affect the membrane dynamics. Whilst stretching <cit.> or compressing <cit.> of the membrane produces wrinkling patterns <cit.>, it has been shown that slowly deforming the membrane, by means of injecting additional fluid into the supporting layer, leads to a dynamics where the fluid pressure is solely balanced by elastic bending forces <cit.>. As the out-of-plane deflection increases, a change in the physical mechanism that dictates the dynamics occurs as elastic stretching becomes the dominant driving mechanism until the system reaches a critical size for which gravity starts to dominate the dynamics. For small systems, an analog to such a predicted transition is observed for a thin perturbed elastic plate resting on a nanoscopic fluid layer, where the restoring elastic bending force is opposed by van der Waals forces leading to an elastohydrodynamic touchdown <cit.> similar to capillary film dewetting <cit.>.
To describe the dynamics theoretically, one can solve the full Navier-Stokes equation in the fluid phase using dynamic boundary conditions at the elastic interface given by the solution of the Föppl-von-Kármán equation <cit.>, using e.g. the immersed-boundary method <cit.>. Viscous flow in thin films can be described by the lubrication theory <cit.> that has been widely used to study different elastohydrodynamic flow phenomena <cit.>. However, not much is known about how elastohydrodynamic flows are affected by the ratio between the geometric parameters that characterize the system as it undergoes large changes while the driving force remain the same.

For instance, when an elastic sheet deforms onto a wall pre-wetted by a thin viscous film, the dynamics of the advancing front are dictated by the local curvature of the interface <cit.>. This elastohydrodynamic relaxation is reminiscent of the capillary spreading of a viscous drop onto a solid substrate <cit.>. Similar to capillary flows, elastohydrodynamic relaxation processes are not only limited to very thin pre-wetted films. In fact, an elastic sheet with zero spontaneous curvature but with an initial shape of a bump (Fig. <ref>) with a height much larger than the pre-wetted viscous film will relax towards a flat equilibrium state. Inevitably, the system must then crossover from a situation where the bump height is larger than the pre-wetted film height, to a situation where instead the pre-wetted film becomes thicker than the bump. Here we investigate how the elastohydrodynamic leveling changes with the ratio between the bump height and the pre-wetted film thickness. In particular, are there different asymptotic regimes, and how do the system transition from one to another? At the nanoscale, thermal fluctuations are expected to contribute and may even dominate the dynamics <cit.>, which we quantify in the leveling dynamics. To answer these questions, we combine numerical solutions of a mathematical model based on the lubrication theory <cit.> with scaling analysis and experiments.



§ MATHEMATICAL MODEL AND NUMERICAL PROCEDURE


We consider the system depicted in Fig. <ref>, where we focus on a system where any influence of gravity can be neglected, i.e., the bump height is smaller than the elasto-gravity length <cit.>. Only situations where the bump height h(x,t) is small compared with its horizontal extent and where the film slopes are small, i.e., ∂ĥ(x,t)/∂ x≪ 1, are considered. We describe the viscous flow between the plate and the solid substrate using lubrication theory <cit.>. When the initial deflection h_i of the elastic plate is small compared to its thickness d, we can neglect stretching and the pressure reduces to p(x,t) = B∂ ^4 ĥ(x,t)/∂ x^4, where B=Ed^3/[12(1-ν^2)] is the bending rigidity of the plate, E is the Young modulus and ν is the Poisson's ratio <cit.>. In addition, the system is a spatially unconfined elastic sheet with the two lateral boundaries being free to move relative to the underlying fluid. Thus, the in-plane compression is suppressed, and bending stresses dominate the relaxation process regardless of the ratio d/h_i. By assuming incompressible flow and imposing no-slip conditions at the two solid substrates, and considering a one-dimensional geometry as there are no variations along the y-direction, one obtains the governing equation for the evolution of the height profile (see e.g. Ref. <cit.>)

    ∂ĥ(x,t)/∂ t = ∂/∂ x[B/12μĥ^3(x,t)∂^5 ĥ(x,t)/∂ x^5 + Γĥ^3/2(x,t)η(x,t)] ,

where μ is the fluid's dynamic viscosity. At small scales, thermal fluctuations can also influence the dynamics, which is described by the last term of Eq. (<ref>). This term mimics the stress generated by thermal fluctuations, originates from an additional symmetric random stress term in the Navier-Stokes equations and is obtained by an integration in the z-direction (for details see <cit.>). The noise term η(x,t), is multiplied by a prefactor Γ=√(k_BT_A/(6μ b)) where k_B is the
Boltzmann constant, T_A is the ambient temperature, b is the width of the plate along the y-direction, and η(x,t) is modelled as a spatiotemporal Gaussian white noise such that ⟨η(x,t)⟩=0 and ⟨η(x,t)η(x',t')⟩=δ(x-x')δ(t-t'), where the ⟨ ⟩ symbols indicate average quantities. We non-dimensionalize Eq. (<ref>) by: X=x/R_i, Ĥ(X,T)=ĥ(x,t)/h_i, T= t Bh_i^ 3 /(12μ R_i^ 6), and Θ(X,T)=η(x,t)√(12μ R_i^ 7/(Bh_i^ 3)). When Γ=0, this non-dimensionalization procedure gives us a parameter-free partial differential equation for Ĥ(X,T). When Γ >  0, the non-dimensional
number N=√(2k_BT_AR_i^ 3/(Bh_i^ 2 b)) appears as a pre-factor in front of the stochastic term, and N^2 measures the ratio between thermal and bending energies. For the macroscopic system provided in our experiment and described in detail below, i.e., T_A = 300 K, h_i=2.5 μm, R_i=20 μm, μ=10^4 Pa s and B=1.3· 10^-12 Nm we get the noise prefactor Γ = 2.5· 10^-13 ms^-1/2 and the energy ratio N=1.75· 10^-6 which is well within the elastic bending dominated regime. However, a transition from a dominant elastohydrodynamic leveling to a dominant stochastic leveling would occur for a system with temperature T_A=300K, membrane perturbation height h_i=10nm and radius R_i=5 μm for a bending modulus B in the range of 10-100 k_BT_A where k_BT_A=4× 10^-21 Nm which corresponds to N in the range of 0-8 <cit.>.

We solve the dimensionless version of Eq. (<ref>) numerically by using a finite element method, and we split it into three coupled equations for: the bump profile H(X,T)=Ĥ(X,T)-ϵ/h_i, the linearized curvature ∂^2 H(X,T)/∂ X^2, and the bending pressure ∂^4 H(X,T)/∂ X^4. These fields are discretized with linear elements and solved by using Newton's method from the FEniCS library <cit.>. For the deterministic case N=0, an adaptive time stepping routine has been used with an upper time step limit of Δ T=0.001 and a discretization in space Δ X ∈ [0.001; 0.01]. For the stochastic case N> 0, we have used a constant time step Δ T=0.001, together with a discretization in space Δ X = 0.0025. At T=0 we impose the initial condition: H(X,T=0) = 1 - tanh(50X^2). We further impose the following boundary conditions at the boundary ∂Ω of the numerical domain: H(X∈∂Ω, T)=H(X∈∂Ω, 0), ∂ ^2 H(X∈∂Ω, T)/∂ X^2=0, and ∂ ^4 H(X∈∂Ω, T)/∂ X^4=0.
 The noise Θ(X,T) is introduced independently at each discrete position and time step using the “random” class with the “randn” Gaussian subclass from the Numpy library <cit.>, with zero mean and a variance 1/(Δ X Δ T). We avoid negative values of Ĥ(X,T) (that might occur in the stochastic case due to the fluctuations), by imposing that when Ĥ(X,T) < 10^-6, it is put back to 10^-6 as in <cit.>.
 To verify the predictions of Eq.(1), we construct an experimental setup which is described in the following section.



§ EXPERIMENTAL PROCEDURE


The experimental setup is composed of a fiber of polystyrene (PS) with a glass-transition temperature T_g, PS≈ 100^∘C deposited on a film of the same polymer supported on a silicon (Si) substrate. These samples are capped by a thin sheet of polysulfone (PSU) with T_g, PSU≈ 180 ^∘C. Sample preparation is carried out as follows: PS fibers (with number-averaged molecular weight M_n = 15.8 kg/mol, and polydispersity index PDI = 1.05, Polymer Source Inc., Canada) are pulled from the melt at 175^∘C using a glass rod. Thin PS films are spin cast from a toluene solution onto 10×10 mm^2 Si substrates, leading to a thickness of 25 to 380 nm, measured using ellipsometry (Accurion, EP3). The films are annealed at 110^∘C for at least 12 hours in vacuum to remove residual solvent and relax residual stresses. The PS fibers are then transferred onto the PS films and the ensemble is heated briefly above T_g, PS. The heating allows the PS to flow, thereby resulting in a bump. Thin PSU films (M_n≈ 22 kg/mol, Sigma-Aldrich) are prepared by spin casting from a cyclohexanone solution onto freshly cleaved mica substrates (Ted Pella, USA). The PSU films have a thickness of ≈ 160 nm, measured using ellipsometry, and are annealed in vacuum at 200^∘C for at least 12 hours. The PSU films are floated onto water and then transferred onto a supporting apparatus (described previously <cit.>), held only by the film edges. These freestanding films can be relaxed to an unstrained state ensuring no in-plane tension. The PSU films are finally transferred onto the PS sample. The part of PSU film at the edges of the Si wafer was then removed using a scalpel blade prior to annealing. This was done to ensure slippage at the boundary between the PSU film and liquid PS layer, thus rendering the relaxation bending dominated as discussed above.

After preparation, the samples were annealed on a hotstage (Linkam, UK) at 130^∘C, which is above T_g, PS but below T_g, PSU. Hence, the PS becomes a viscous liquid while the capping PSU film remains an elastic solid, thus realizing the system illustrated in Fig. <ref>. The height profile is imaged during annealing using optical microscopy with a red laser line filter (λ = 632.8 nm, Newport, USA), which creates interference fringes in the region of the bump, as shown in Fig. <ref>(a), due to the light that is reflected from the Si substrate. It is clear from these fringes that the initial fiber and resulting flow are one-dimensional over length scales that are many times the width of the perturbation itself. Each interference fringe corresponds to a change in height of λ/(4n), where n ≈ 1.57 is the average index of refraction of the two polymers that make up the sample (n_PS = 1.53 and n_PSU =1.61). This allows the bump profile h(x,t) to be reconstructed by fitting a polynomial to the fringe data, as shown in Fig. <ref>(b). Such profiles can then be used to track the leveling dynamics, and to extract in particular the evolution of the height h_0(t) of the bump with time for various initial geometries.



§ RESULTS



 §.§ Elastohydrodynamic leveling




We first start by investigating the elastohydrodynamic leveling in the absence of thermal fluctuations (N=0). In Fig. <ref> we show the numerical solutions of the dimensionless version of Eq. (<ref>) for N=0 and we can see that the aspect ratio h_i/ϵ controls both the time scale for leveling and the detailed features of height profile. The smaller h_i/ϵ, the faster the dimensionless leveling process. Also, the dip created near the advancing front of the perturbation is enhanced both in magnitude and lateral extent for smaller h_i/ϵ. We remark that for each initial aspect ratio there is a transition period of a few numerical time steps preceding the onset of the leveling process. This part of the data is not included in Fig. <ref> as it is considered to depend on the initial condition, but does not influence the later dynamics.  For h_i/ϵ=43, the numerical height profiles are further compared with our experiments, which are found to be in good agreement. We recall here that the elastic plate is floating on the liquid film and has edges that are free to move. Therefore, the pressure contribution from bending still largely dominates any contribution from stretching and Eq. (<ref>) is still valid.


We now turn to a scaling analysis of Eq. (<ref>) for N=0. When h_0(t)/ϵ≪ 1, the equation can be linearized and reduces to 12μ∂ h/∂ t=Bϵ^3∂^6h/∂ x^6 and we deduce the long-term scaling for the temporal evolution of the horizontal length of the bump: R(t)∼ [Bϵ^3t/(12μ)]^1/6. Since there is area conservation in the (x,z)-plane, we assume R(t)h_0(t) to be constant, that is evaluated to R_ih_i at t=0. By combining these scaling relations we get for h_0(t)/ϵ≪ 1

    h_0(t)/ϵ∼(τ/t)^1/6 ,

where τ = 12μ h_i^ 6R_i^ 6/(B ϵ^9) is the characteristic time scale for the bending-driven leveling dynamics. As we operate within the regime where bending dominates over stretching, a similar result is obtained by considering the force balance between the viscous and bending forces <cit.>. Also if we include isotropic stretching due to clamped boundaries a similar scaling law appear, but now with an additional logarithmic term, R(t)∼ (t/log(t))^1/6 <cit.>. However, when h_0(t)/ϵ≫ 1 we must match the curvature of a traveling-wave solution localized near the advancing front with the quasi-static solution to obtain the correct scaling <cit.>, i.e., constant pressure in the bump, leading to <cit.>

    h_0(t)/ϵ∼(τ/t)^2/17 .

By balancing the two asymptotic predictions above, we expect the crossover between them to occur around t/τ≈1. In addition, these asymptotic regimes suggest that h_0(t)/ϵ is essentially a function of t/τ only, independent of the value of h_i/ϵ in particular.

In order to test our scaling predictions, we compute numerical solutions of the dimensionless version of Eq. (<ref>) for N=0, with h_i/ϵ∈[10^-2, 10^3] and extract h_0(t)/ϵ as a function of t/τ. These numerical results are plotted in Fig. <ref> and compared to the experimental data. For each sample, the experimental data is matched to the numerical data through one fitting parameter α in front of the time scale τ. The values of the fluid viscosity and elastic Young modulus are highly sensitive to the temperature in the experiments, and we estimate them to be μ≈10^4 Pa s <cit.> and E≈2.6 GPa <cit.>, respectively. Since all experiments were carried out at the same temperature and with the same polymer, sample-to-sample variations in τ result only from uncertainties in the geometrical parameters h_i, R_i, d, and ϵ. The obtained α values are 0.13, 0.7, and 1.3 for the three samples and each of these values are reasonably close to unity. More importantly, the sample-to-sample variations in α do not exceed a factor of 10, which is well within the expected relative error arising from the high sensitivity of τ to the geometrical parameters. The general agreement between the experimental data and the numerical predictions is good, over about 5 orders of magnitude in t/τ. The systematic early-time tail in the experimental data might be attributed to the initial compressive thermal stresses in the elastic layer, which arise due to the rapid heating of the samples from room temperature to T = 130^∘C, which relax prior to leveling and the time needed for the initial shape to enter the asymptotic regime.

The master curve in Fig. <ref> confirms that h_0(t)/ϵ is a function of t/τ. Furthermore, the two scaling regimes predicted above are indeed present, with pre-factors close to unity, and the crossover between the two being located near t/τ≈1 as predicted. Any bump that initially starts in a thin pre-wetted film regime h_0(t)/ϵ≫1 will eventually crossover to a thick-film regime h_0(t)/ϵ≪ 1, with the corresponding power laws in time. As a final remark, a similar combination (not included here) of numerical simulations and scaling analysis can been performed for an axisymmetric geometry, leading to h_0(t)∼ t^-2/11 for h_i/ϵ≫ 1, and h_0(t)∼ t^-1/3 for h_i/ϵ≪ 1.



 §.§ Stochastic leveling

Next we investigate the leveling process when it is dominated by thermal fluctuations (N> 0). As shown in Fig. <ref>, the numerical solutions suggest that the aspect ratio h_i/ϵ is again essential, as it sets the time scale for leveling where the smaller h_i/ϵ, the faster the dimensionless leveling process. Moreover, by comparison with the deterministic (N=0) case in Fig. <ref>, the stochastic (N> 0) profiles exhibit spatiotemporal fluctuations and adopt different average shapes and leveling dynamics.



To go further, we propose a scaling analysis of Eq. (<ref>), inspired by Ref. <cit.>. We consider specifically the N≫1 limit, for which the thermal fluctuations are the dominant driving contribution to the dynamics and we assume that we can neglect the bending term so that Eq. (<ref>) reduces to ∂ h/∂ t= Γ∂[(ϵ+h)^3/2η]/∂ x. We consider the average quantities ⟨ h_0(t)⟩ and ⟨ R(t)⟩, where we invoke the ∼ (tx)^-1/2 scaling <cit.> for the root mean square value of the averaged noise over a space interval x and a time interval t. By assuming that the average area conservation in the (x,z)-plan can be expressed as ⟨ h_0(t)⟩⟨ R(t)⟩∼ h_iR_i, we get

    ⟨ h_0(t)⟩/ϵ[1+⟨ h_0(t)⟩/ϵ]^3∼τ_Γ/t

where τ_Γ = 6μ h_i^ 3R_i^ 3 b/(k_B T_Aϵ^4) is the characteristic time scale for the stochastic leveling dynamics. Interestingly, Eq. (<ref>) describes a complete crossover between two asymptotic regimes in the stochastic leveling dynamics: for ⟨ h_0(t)⟩/ϵ≫ 1, we obtain ⟨ h_0(t)⟩/ϵ∼(τ_Γ/t)^1/4, and thus we recover ⟨ h_0(t)⟩∼ t^-1/4 <cit.>; while for ⟨ h_0(t)⟩/ϵ≪ 1, we get ⟨ h_0(t)⟩/ϵ∼τ_Γ/t. We expect the crossover between the two asymptotic regimes to occur around ⟨ h_0(t)⟩/ϵ≈ 1, i.e., around t/τ≈1/8.

In order to test the prediction in Eq. (4), we compute the numerical solution of the dimensionless version of Eq. (<ref>) for 5≤ N≤ 8, with h_i/ϵ∈[10^-1, 10^2]. By averaging over a minimum of 30 realizations, we can extract ⟨ h_0(t)⟩/ϵ as a function of t/τ_Γ and the results are plotted in Fig. <ref>. The data from the numerical solutions is in good agreement with Eq. (<ref>) for all ⟨ h_0(t)⟩/ϵ, with no adjustable parameter. Our results highlight that Eq. (<ref>) gives an accurate prediction of the stochastic leveling dynamics and show that the missing pre-factor is close to unity. Finally, in order to further highlight the underlying self-similarity associated with each of the two asymptotic regimes, the insets of Fig. <ref> show the corresponding bump height profiles rescaled according to  Eq. (<ref>). In each asymptotic regime the height profiles collapse onto a universal shape which confirms the overall self-similarity in the leveling dynamics.



§ CONCLUSION

We have described the elastohydrodynamic and stochastic leveling of an elastic plate placed atop a viscous film. By combining numerical solutions, scaling analysis, and experiments, we identified various canonical regimes. Our results highlight the importance of the driving mechanism, either by elastic bending of the plate or thermal fluctuations, and the influence of the aspect ratio of the bump height to the film height. For each of these two driving mechanisms, a crossover between two distinct asymptotic regimes is controlled by the aspect ratio. These findings can be helpful to explain elastohydrodynamic leveling dynamics found in biological, engineering, or geological processes.


The financial support by the Research Council of Norway (project grant
263056), the Natural Science and Engineering Research Council of Canada and the Joliot Chair of ESPCI Paris is gratefully acknowledged. The authors also thank Tak Shing Chan, Miroslav Kuchta, Corentin Mailliet, René Ledesma-Alonso, Maxence Arutkin, Elie Raphaël, Howard Stone, Pauline Olive, Katelyn Dixon, Paul Fowler, Mark Ilton and David Dean, for interesting discussions and preliminary works in different geometries.


unsrt

