Viscous Fluid Instabilities Under An Elastic Sheet by Maria Khomenko A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) May 2010 c Maria Khomenko 2010 Abstract This thesis considers the flow of thin fluid film between an elastic sheet and a rigid plane. We derive a mathematical model for the flow from the NavierStokes equations using the lubrication approximation and develop numerical and similarity solutions to this problem. An experimental apparatus was developed to investigate this phenomenon, and the results of the mathematical model were compared with experimental data. Chapter 3 examines the evolution of a fixed fluid volume under gravitational forces on a horizontal plane. The evolution of the fluid mass profile and the progression of the fluid front are determined from the numerical solutions, as well as experimentally. The favourable comparison between the numerical solutions and the experimental results establishes the validity of the model. Chapters 4-5 considers the evolution of a thin fluid flow under an elastic on an inclined plane. We establish a traveling wave solution for this flow. A linear stability analysis yields the criterion for the existence of unstable modes and establishes the growth rate and wavelength of the most unstable mode. Instability is promoted by increasing the inclination of the plane. For low angles, the numerical and experimental growth rates were in good agreement, while the wavelengths were experimentally of the same order and numerically computed wavelengths had little variation. The long term behaviour of the fluid front is studied analytically via a similarity solution in Chapter 6. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication xi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Derivation of the Governing Equations . . . . . . . . . . . 2.1 Equation for Velocity . . . . . . . . . . . . . . . . . . . . . 2.2 Imposing the No-slip Boundary Condition . . . . . . . . . . 4 4 6 3 Horizontal Flow (Special Case with α = 0) . . . . . . . . . . . . . . . . 3.1 Non-dimensional Variables . . . . . . . . . . . . . . 3.2 Numerical Analysis . . . . . . . . . . . . . . . . . . 3.2.1 Initial Condition . . . . . . . . . . . . . . . . 3.2.2 Steady Initial Shape . . . . . . . . . . . . . . 3.2.3 Numerical Solution of the PDE . . . . . . . . 3.3 Experimental Results . . . . . . . . . . . . . . . . . 3.3.1 Experimental Effects of the Bending Stiffness 3.3.2 Numerical-Experimental Comparison . . . . . 3.4 Similarity Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 8 10 11 12 14 17 19 25 29 4 Inclined Plane: Numerical Model . . . . . . . . . . . . . . . 4.1 Non-dimensional Form . . . . . . . . . . . . . . . . . . . . . 35 35 iii Table of Contents 4.2 4.3 The Propagating Wave Front . . . . . . . . . . . . . . . . . Time Dependent Numerical Analysis . . . . . . . . . . . . . 4.3.1 Solution with Constant Flux Boundary Conditions . 4.3.2 Stability Analysis . . . . . . . . . . . . . . . . . . . 4.3.3 Stability Analysis: Small Wavenumber Limit q 1 . 4.3.4 Solution with Constant Volume Boundary Conditions 5 The Elastic Experiment . . . . . . 5.1 Dimensional Variables . . . . . . 5.2 Contour Data Processing . . . . 5.2.1 Growth Rate . . . . . . . 5.2.2 Wavelength . . . . . . . . 5.2.3 Velocity . . . . . . . . . . 5.3 Alternative Dimensional Scaling 5.4 Discussion of Experimental Errors . . . . . . . . 58 58 61 62 65 66 68 70 . . . . . . . . . . . . . . . . . . . . . . . 72 7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . 79 79 80 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6 Similarity Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 41 41 43 48 54 Appendices A Numerical Details . . . . . . . . . . . . . . . A.1 Horizontal Flow: Initial Condition . . . . A.2 Numerical Solution of the PDE . . . . . . A.2.1 Time Discretization . . . . . . . . A.2.2 Space Discretization . . . . . . . . A.2.3 Numerical Boundary Conditions . A.2.4 Extending Numerics for Solving the A.3 Inclined Plane: Solution of the ODE 4.11 A.3.1 Growth Rate Eigenvalue Solution . . . . . . . . . . . . . . . . . . . PDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 84 85 85 86 88 89 90 91 B Experimental Set Up . . . . . . . . . . . . . . . . . . . . . . . 92 C MATLAB Code . . . . . . . . . . . . . . . . . . . . . . . . . . 94 iv List of Tables 5.1 5.2 5.3 5.4 5.5 The values for hc = 3mm . . . . . . . . . . . . . . . . . . . . Growth rates . . . . . . . . . . . . . . . . . . . . . . . . . . Characteristic time . . . . . . . . . . . . . . . . . . . . . . . Characteristic length with numerical wavenumber and wavelength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The values for hc = 3mm . . . . . . . . . . . . . . . . . . . . 59 64 64 66 69 v List of Figures 2.1 3.1 3.2 3.3 Approximate shape of a two dimensional cross section of the fluid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fluid on a horizontal surface . . . . . . . . . . . . . . . . . . Pre-wetted front . . . . . . . . . . . . . . . . . . . . . . . . . Cross-sectional view. The fluid is between two rollers, which have radius R and are located distance d from the center. The fluid follows the circular shape of the roller until a ‘contact point’ at angle θ. . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Solutions for different θ’s (as labeled in the legend, in radians). 3.5 Numerical simulation of a horizontal fluid bump spread. . . . 3.6 (a) The height decrease with respect to time. (b) Front propagation with respect to time. . . . . . . . . . . . . . . . . . 3.7 The raw experimental result from the first frame. The shape symmetrized for purposes of numerical analysis. . . . . . . . 3.8 The tallest peak corresponds to the fluid profile 2 seconds after the release and the last profile is three minutes after the release. Three profiles in between are included. . . . . . . . . 3.9 The fluid shape for several experimental runs. The discrepancy in height is due to different volumes. . . . . . . . . . . 3.10 The initial experimental shape of the fluid (dotted line) approximated by the solutions of the tension equation (dashed) and the bending stiffness equation (solid line). . . . . . . . . 3.11 Comparison of numerical simulations starting with the numerical initial conditions governed by bending stiffness (solid line), and by tension (dashed line) and non-dimensional symmetrized experimental shape (dotted line) after 4 minutes. . 3.12 Comparison of numerical simulation (solid line) and experimental results (dotted line) at different times. Slight discrepancies in the volume of the fluid are due to experimental errors when measuring the fluid height. . . . . . . . . . . . . 5 9 10 13 14 15 16 18 19 21 22 24 26 vi List of Figures 3.13 Comparison of numerical simulations with bending stiffness (solid line) and with tension (dashed line) at time 9 seconds. Both use experimental steady state shape as initial conditions. 3.14 Top: Numerical, with experimental initial conditions, simulation (solid line) and experimental (dotted line) front positions with respect to time on a logarithmic scale. Bottom: Numerical, with experimental initial conditions, simulation (solid line) and experimental (dotted line) fluid heights with respect to time on a logarithmic scale. . . . . . . . . . . . . 3.15 Position of the front as a function of time. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. 3.16 Height of the fluid as a function of time. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. 3.17 Height of the fluid as a function of the position of fluid the front. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. . . . . . . . . . . . . . . . . . . . . . 4.1 4.2 4.3 4.4 4.5 4.6 Top: Different δ’s. Bottom: Different D values. . . . . . . . Top: Different δ values. Bottom: Different D values. . . . . Top: Leading order solution used as h0 . Bottom: Solution of the eigenvalue σ as a function of q. The mode that will grow preferentially is the q corresponding to the maximum value of σ(q). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Growth rate dependence on δ value. . . . . . . . . . . . . . . A comparison of the small order q approximation with the full solution for growth rate, σ(q), at inclination of α = 60o (D=0.1 and δ = 0.05). . . . . . . . . . . . . . . . . . . . . . The approximated growth rate for small q values for different values of D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 28 32 33 34 40 42 47 48 52 53 vii List of Figures 4.7 4.8 4.9 5.1 5.2 5.3 5.4 5.5 5.6 6.1 6.2 Numerical solution of the PDE (4.10) with the constant volume boundary conditions. . . . . . . . . . . . . . . . . . . . Tracking the shape of the wave front over a long time (up to t=50). This tracks the height and the position of the peak with time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Position of the wave front against time (dimensionless form). Top: Numerical solution (δ = 0.05) with the indicated height used as hc = 0.5853. Calculated for time = 2 which is 9.8 seconds when converted to dimensional units for this image, but this conversion depends on tc and may vary for different inclinations. Bottom: Experimental behaviour with corresponding hc = 2.987 at 10 seconds after the release. This time may vary for individual experiments from 10 to 15 seconds . . . . The comparison of experimental and numerical results with angle of inclination of 60o . With tc = 4.9s and hc = 3mm. Top: Comparison 10 seconds after the release with constant volume numerical solution. Fluid is gradually transforming into the shape of the traveling wave front. Bottom: Comparison 1 minute after the release with the constant flux numerical solution. Note, the origin of the horizontal scale has been shifted. . . . . . . . . . . . . . . . . . . . . . . . . . . . The contours of fluid fronts at inclination of 78o from time 0 to 133 seconds. . . . . . . . . . . . . . . . . . . . . . . . . . Measuring growth rates . . . . . . . . . . . . . . . . . . . . . Growth rates . . . . . . . . . . . . . . . . . . . . . . . . . . Top: The distance traveled by the wave front from some initial level (the initial level can be chosen arbitrarily, since we are interested in the slope). The black solid line indicates the shape of the fitted polynomial. Bottom: The speed of the wave front, as calculated from the derivative of the distance traveled. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Position of the front with respect to time. Experimental results (dashed line), numerical simulation (solid line) and similarity solution from equation 6.7 (dotted line). . . . . . . . . Fluid height at the front with respect to time. The similarity solution is given by (6.8). . . . . . . . . . . . . . . . . . . . . 55 56 57 60 61 62 63 65 68 75 77 viii List of Figures 6.3 Fluid height at the front with respect to front position. The similarity solution is given by (6.9). . . . . . . . . . . . . . . 78 A.1 The simulation of the horizontal flow with three δ values varying from 0.01 to 0.1. Each simulation run is graphed at four different times. . . . . . . . . . . . . . . . . . . . . . . . . . 89 B.1 Experimental set up . . . . . . . . . . . . . . . . . . . . . . 93 ix Acknowledgements I would like to express my deepest gratitude to my co-supervisors, Neil Balmforth and Anthony Peirce, for their expert guidance, patience and unwavering support. This work would not have been possible without their understanding and encouragement. I am also grateful to the Department of Mathematics, the Institute of Applied Mathematics for helping and the Complex Fluids Lab for providing research facilities as well as a welcoming environment. I would like to thank my undergraduate professors, James Forrest, Zoran Mišković and Francis Poulin, who took the time to work with me and helped me to be where I am today. Finally, I want to thank Karsten Chipeniuk, who was always there for me. Funding for this project was provided by NSERC CGS fellowship. x ”Those who mind don’t matter and those who matter don’t mind.” ∼ Dr. Seuss ∼ To those who matter... xi Chapter 1 Introduction Studies of thin fluid film behaviour date back to the 1800s, beginning with the works of Young and Laplace [3]. Today, thin fluid flows have a number of applications in various fields. These include, but are certainly not limited to, the study of spin-coating processes and microchip production in the electrical engineering industry, the lining of mammalian lungs in biology, and formation of water rivulets along the wing of an aircraft [1, 8]. Traditional analysis of these problems assumes that the top surface layer of the film is free, and that the shape of the fluid is governed by the surface tension while the motion is driven by gravity. This leads to the governing equation being a fourth order partial differential equation (PDE), which depends on the surface tension γ of the fluid and the gravitational acceleration, g = 9.81m/s−2 . There are numerous works, both theoretical and experimental, which analyze this surface tension model [1, 7]. The work in this thesis takes an unconventional approach to the classical model of thin fluid flow. The top layer of the fluid is no longer treated as a free surface governed by the surface tension. Instead it is now covered with an elastic membrane (or plate). This effectively changes the top boundary condition to a no-slip condition and the shape of the fluid is now governed by the bending stiffness of the membrane. Hence, the governing equation gains two extra derivatives, making the model a sixth order PDE. While the two models share a number of common attributes, the physical problems have qualitative and quantitative differences, as well as different families of practical applications. Behaviour of a fluid under a thin film is of great interest in engineering industries, such as working with silicon wafers [11] and semiconductor manufacturing [13]. Moreover, thin films forming on the surface of a liquid have a number of geological applications, such as the motion of the molten interior of a lava flow under the solidifying crust [13]. Introducing the thin elastic film, and effectively removing the free surface can have a significant effect on the flow. The fluid cannot ignore the presence of the new boundary condition and follow the motion dictated by the surface tension, so it must adjust its flow pattern to balance its natural 1 Chapter 1. Introduction behaviour with the interference from the film. Some of the known effects are formation of wrinkles [6] or buckling patterns [13]. These phenomena are not necessarily strong enough to alter the overall fluid behaviour, and yet the effects are of sufficient magnitude to be of interest in the applications mentioned above. The mathematical modeling of a fluid flowing under constraints can be extended further still. There have been a number of studies in the past few decades modeling the flow through an elastic solid, such as magma uptake through the earth during a volcanic eruption [12], and industrial application of enhanced oil recovery. In the latter, a fluid is used in the formation of hydraulic fracture [2]. This formulation introduces different boundary conditions. However, the general governing equations are still based on lubrication theory and the general shape of the steady state solution is similar. We will look at how the extra stiffness of a thin film affects the flow of a fluid down an inclined plane. The qualitative analysis is similar to previous works dealing with surface tension; we still model the problem though lubrication theory and use a thin flow approximation to arrive at the governing equation. However, while the two problems share a number of similar qualitative features, the quantitative results differ by orders of magnitude. The numerical computations are also compared with experimental results. Particularly, we use light intensity to measure the height of the fluid and compare its general shape to that predicted by the numerical model. We also monitor the progression of the fluid down the inclined plane, which allows us to determine growth rate and wavelength of instabilities, which form during the flow, as well as the velocity of the flow. The long term behaviour of the fluid is also studied though a similarity solution. This is a commonly used approach, which is an asymptotic reduction of the governing equations to a simpler form, which is often a single non-linear PDE [11]. The study is presented in two parts, each containing a theoretical and experimental component. The first study looks at the flow on a horizontal plane, while the second introduces inclination to the problem. Setting the angle of inclination to 0 slightly simplifies the numerical aspects of the modeling and significantly increases our control of the experiment. The results are more straight forward to compare and the measurements are more accurate. Thus, the first part establishes the validity of the model, while the 2 Chapter 1. Introduction second looks at more interesting behaviour of the fluid, such as instability evolution. The experimental component measures the actual physical behaviour of a viscous fluid covered by an elastic sheet. The fluid is allowed to move either on a horizontal surface or down an inclined plane and its shape at subsequent time steps is observed. When there is inclination present, the fluid will tend to form instabilities (or fingers) as it is flowing down a plane. We are interested in observing the effects of introducing the elastic sheet, and in particular how the behaviour of the fluid differs from when it is allowed to move under surface tension. Data is collected at several different angles of inclination and analyzed. The next section, Chapter 2 introduces the general form of the model used throughout the study and derives the governing thin fluid film equation. Chapter 3 contains all the analysis (theoretical and experimental) pertaining to the flow on a horizontal surface. The comparison between experiment and theory in this section justifies the validity of the model. Chapter 4 applies the model to the flow on an inclined plane and examines the theoretical behaviour of the instabilities. Lastly, Chapters 5 and 6 compare the numerical predictions from Chapter 4 with the experimental data and look at the long term behaviour of the fluid. 3 Chapter 2 Derivation of the Governing Equations 2.1 Equation for Velocity Consider the problem of a viscous fluid flowing between an elastic sheet and an inclined plane. As the fluid moves along the surface, it spreads out and starts behaving as a thin fluid film with boundary conditions which depend on the bending stiffness of the elastic sheet. At this point, thin film flow approximations become applicable [7]. When dealing with the modeling of fluid flow, the starting point is the Navier-Stokes equation: 1 µ ∂~u + (~u · ∇)~u = − ∇P + ∇2~u + g sin αî − g cos αk̂, ∂t ρ ρ (2.1) where ~u = (u, v, w) = (~v , w) is the fluid velocity, ρ is the fluid density, P is the pressure, and µ is the viscosity, which is a measure of frictional forces present in the fluid. The last two terms are the downhill and the normal components of gravitational acceleration, g, respectively. The fluid film considered is also assumed to be incompressible, which can be expressed as the divergence-free condition, ∇ · ~u = 0. (2.2) The shape of the fluid is shown in Figure 2.1. The features of interest are the capillary ridge immediately behind the fluid front, the contact point of the fluid front with the surface and the fluid thinning behind the front. These matters are addressed in Chapter 3. 4 2.1. Equation for Velocity Elastic sheet Fluid film v h(x) Figure 2.1: Approximate shape of a two dimensional cross section of the fluid. To model the physical problem we assume that the fluid is thin, that is, the height of the fluid, h(x, y), is much smaller than the horizontal scale. Thus, for a thin film fluid flow, we can define a ratio = h/L, where L is the horizontal scaling. Then the horizontal derivatives of the velocity with respect to x, i.e. tangent to the plane, are much smaller than the vertical derivatives with respect to z (normal to the plane). That is, in non-dimensional form, ∂ 2 |~v | ∂ 2 |~v | . ∂x2 ∂z 2 Therefore, the in-plane derivatives of the velocity are negligible in the momentum equation. Also, since the flow has a very low Reynolds number (∼ 0.1 given experimental parameters), the inertial terms are small in comparison to the viscous forces and equation (2.1) reduces to 2 ∂ 2~v + ρg sin αî − ρg cos αk̂. (2.3) ∂z 2 In vector form, separating out the in-plane components and the normal component, 0 = −∇P + µ ∂ 2~v + ρg sin αî, ∂z 2 ∂P = −ρg cos α, ∂z ∇P = µ (2.4) (2.5) where equation (2.4) is now in 2 dimensions, so ∇ = (∂x , ∂y ). 5 2.2. Imposing the No-slip Boundary Condition Equation (2.5) can be used to solve for pressure. The upper layer boundary condition at z = h(x, y) for a thin fluid film is approximated to be P (h) = B∇4 h, where B is the bending stiffness of the elastic sheet. The pressure as a function of z is P (z) = −ρg(z − h) cos α + B∇4 h. (2.6) Substituting this result back into equation (2.4) we obtain an expression for the velocity: ~v = (∇P − ρg sin αî) z2 ~ + C, ~ + Az 2µ (2.7) ~ and C ~ are to be determined by the boundary condiwhere the constants A tions. 2.2 Imposing the No-slip Boundary Condition The no-slip boundary condition is present at the surface layer z = h(x, y, t), so the tangential velocity is zero and the normal velocity is the change in height (ht ). The normal to the surface is n̂ and tangents are t̂1 and t̂2 , so that un = n̂ · (u, v, w) = ht , t̂1 · (u, v, w) = 0, t̂2 · (u, v, w) = 0. (2.8) (2.9) The normal is n̂ = (−hx , −hy , 1). However, when accounting for the difference between the component in the z direction and the changes in the x direction, the term hx is of order . Thus, the normal velocity component is dominant here, and in-plane components are negligible. Through similar treatment of the tangential components, the boundary surface conditions are u, v ≈ 0, w = ht . (2.10) 6 2.2. Imposing the No-slip Boundary Condition The boundary adjacent to the plane has a no-slip condition as well and there is no penetration. That is, the volume is conserved and the tangential component vanishes: ~v = 0. ~ = 0. Similarly, it Combining this with equation (2.7), z = 0 implies C ~ = −(∇P − ρg sin αî) h . Equation (2.7) follows that z = h gives the value A 2µ now becomes ~v = (∇P − ρg sin αî) z (z − h). 2µ (2.11) Averaging over the vertical direction removes the z-dependence from ~v , Z 1 h h2 h~v i = ~v dz = − (∇P − ρg sin αî). (2.12) h 0 12µ Since there is no fluid loss, the conservation of mass holds and it requires that ∂h + ∇ · (hh~v i) = 0. (2.13) ∂t The explicit form of ~v and P are already known from equations (2.6) and (2.12). Hence, the conservation of mass can be re-expressed as a PDE for height, h: h i 1 ∂h = ∇ · Bh3 ∇∇4 h + ρgh3 cos α∇h − ρgh3 sin αî . ∂t 12µ (2.14) Equation (2.14) will be referred to as the dimensional thin film equation. The height is the primary unknown of the problem and all the other parameters µ, B, ρ, g and α are known physical quantities. For further numerical analysis, the equation must be converted to its non-dimensional form to ensure the results are as general as possible. The non-dimensional scalings used are dependent on the physical problem. In particular, two different scalings are used, corresponding to the two different physical situations. The first problem describes the thin fluid flow on a horizontal surface with elastic boundary conditions. The second problem describes the thin fluid flow on an inclined plane. 7 Chapter 3 Horizontal Flow (Special Case with α = 0) It is difficult to deal with a high order differential equation, such as equation (2.14). It is more difficult still to implement an experiment depending on many parameters and variables, each of which introduces experimental errors. Hence, in this section, we consider the fluid flow to be perfectly horizontal and so we take the angle of inclination, α, to be zero. This simplified case removes an extra parameter, thus removing an extra dimensional group, and allows us to concentrate on the effect of the bending stiffness on the fluid wavefront. In order to produce a meaningful initial condition, two cylindrical rollers are used to compress the fluid into an elongated bump. The experiment becomes much more manageable with this simplification, as the process of measuring the fluid height becomes much easier and removal of the extra parameter reduces uncertainty. The inclined plane is considered in the next chapter. 3.1 Non-dimensional Variables The full governing equation of the motion (derived in Chapter 2) is h i 1 ∂h = ∇ · Bh3 ∇∇4 h + ρgh3 cos α∇h − ρgh3 sin αî . ∂t 12µ (3.1) Introducing the α = 0 simplification removes the last term and reduces the second term, while leaving the bending stiffness term unaltered. In vector form, the equation becomes ∂h 1 = ∇ · Bh3 ∇∇4 h + ρgh3 ∇h . ∂t 12µ (3.2) 8 3.1. Non-dimensional Variables The next major simplification comes from the symmetry of the problem. Figure 3.1 shows a typical starting fluid bump. Since there is no inclination, when the fluid is released, the wave front continues to move independently of the y position. Effectively this reduces the motion to two dimensions and all the important information is contained in a single x − z cross-section. In turn, the governing equation (3.2) is reduced to the 1-dimensional equation z x 1 3 ∂h = Bh hxxxxx + ρgh3 hx x . (3.3) ∂t 12µ y Introducing the scaling transformation (x̄, h̄, t̄) = (x/xc , h/hc , t/tc ), it is required that the time derivative in equation (3.3) is Figure 3.1: Fluid on a horizonof the same order as the right hand side: tal surface 1 1 1 1 hc = Bh3c hc , tc 12µ xc xc x4c hc 1 1 1 = ρgh3c hc . tc 12µ xc xc This yields the required time and length scalings: xc = B ρg 1/4 , tc = 12µx6c . Bh3c (3.4) The thin film equation can be reduced to its non-dimensional form with these scalings. The bars on the non-dimensional variables are dropped for simplicity, yielding the resulting equation ht = h3 hxxxxx + h3 hx x . (3.5) From (3.4) we observe that there are three unknown characteristic parameters, xc , tc and hc , and only two equations relating them to each other. All other quantities are known physical properties of the fluid. Thus, one of these parameters has to be determined and in order to obtain a comparison between the solution and the experimental results its value is chosen to agree with the experiment. The characteristic time, tc , is a general measure of how fast the fluid moves. The characteristic length, xc , is a general measure of how wide the 9 3.2. Numerical Analysis fluid front spreads, and the characteristic height is the initial height of the fluid film. Of the three characteristic variables, the last one is the only one that is associated with a specific physical value which can be measured. Hence, hc is chosen as the pre-defined parameter and is assigned the corresponding experimental height value (typically an order of magnitude of 5-10 mm). 3.2 Numerical Analysis The problem is modeled by a sixth order non-linear PDE (equation (3.5)) which is degenerate and very stiff, making the problem a non-trivial one to solve. One of the present mathematical difficulties is computing the solution near the contact line. Particularly, Elastic Plate at the point where the edge of the h(x,t) δ fluid makes contact with the surface, the fluid height drops to zero and Figure 3.2: Pre-wetted front the use of a no-slip boundary condition, which requires that the parallel velocity components are zero. Mathematically, this would cause the fluid height at the contact line to go to infinity and make the slope of h infinitely steep, thus, creating a singularity. This classic problem is present in the surface tension model as well and there have been numerous works that have looked into the specifics of the contact line [15, 4]. Numerically, this issue is resolved by assuming a pre-wetted film around the fluid. That is, when the fluid moves, it is peeling the elastic sheet off a thin existing layer of fluid. This extra simplification is also consistent with the experimental setup, where the domain of the flow is in fact pre-wetted. However, it is worth noting that in the analysis of the similarity solution, Section 3.4 and Chapter 6, this assumption may be relaxed. Particularly, because the boundary conditions do not directly affect the outcome. We also relax this condition in the analytic stability analysis, Section 4.3.3, where the analytic derivation does not necessarily require the boundary condition at the contact line. There are two ways to look at the numerical modeling. One is to consider a bump in the centre of the domain. At time 0 the bump is released and 10 3.2. Numerical Analysis is allowed to spread. The boundary conditions on each side are identical, namely the height is kept at a precursor film thickness δ 1 and the derivatives are zero. On a domain [0, L] the boundary conditions are: h(0, t) = δ, h0 (0, t) = 0, h00 (0, t) = 0, h(L, t) = δ, h0 (L, t) = 0, h00 (L, t) = 0. (3.6) Alternatively, one could exploit the symmetry of the spread and only consider half the domain. Then the right hand side conditions would remain the same, but left hand side would correspond to the centre of the bump. At the centre, the height is continuously decreasing and cannot be predefined, but since it is known that the flow is symmetric, all the odd derivatives at that point must be zero. The boundary conditions in this case are: h0 (0, t) = 0, h000 (0, t) = 0, h00000 (0, t) = 0, h(L, t) = δ, h0 (L, t) = 0, h00 (L, t) = 0. (3.7) Working with a smaller domain allows for greater numerical accuracy, while still keeping the computing time reasonable. The two methods do in fact yield the same results, provided they start with the same initial conditions. The details of the numerical methods employed are described in Appendix A. 3.2.1 Initial Condition The motion of the fluid strongly depends on the pressure exerted by the fluid on the elastic sheet. The pressure at a point is in turn highly dependent on the shape of the fluid and how much mass is close enough to contribute to the force at a given point on the surface. Thus, the initial shape of the fluid bump has a very strong effect on how the fluid will move. Flows starting with significantly different boundary conditions will move at different speeds and may have different shapes when compared on a similar time scale. This makes the numerical problem very sensitive to the choice of initial conditions. One important factor to consider is that the fluid volume 11 3.2. Numerical Analysis of the simulation and the experiment must agree exactly in order for the two to be compared. Thus, in order to compare the numerical and experimental results, the initial conditions of the simulation must be calculated based on the experimental parameters (such as fluid volume, initial height and distances). This is an involved calculation, which depends on multiple parameters and can become rather cumbersome. A more effective approach is then to take the initial conditions directly from the initial shape of the fluid immediately following the release. This guarantees that the simulation and the experiment are tracking the identical motion. The comparison is included in Section 3.3. A purely computational solution (with initial conditions determined numerically) will not necessarily correspond to the same physical initial condition as the experiment and we cannot expect perfect agreement from a rigorous comparison of the two. Yet, they will still share the same attributes and it is still worthwhile to analyze the numerical model. Consistency between the general behaviour of the numerical model and the experimental results would also solidify the comparison analogy and discrepancies may suggest possible experimental errors. 3.2.2 Steady Initial Shape For calculating numerical initial conditions consider the equation P ρg h+ . (3.8) B B This equation dictates the shape of the fluid subject to boundary conditions from the rollers, which hold the fluid in place, and symmetry. Symmetry demands that the odd derivatives at the centre (x = 0) are zero. The transition between the segment of the elastic sheet which is in contact with the roller and the segment between the rollers should be smooth. That is, the derivatives are continuous and at the point of contact the derivative of the sheet curvature is equal to the curvature of the roller circle. Define this contact point in terms of an angle θ, which is the angle between the point of contact and the lowest point on the circle (Figure 3.3). hxxxx = At the right hand contact point the derivatives are: dh = − tan θ, dx (3.9) 12 3.2. Numerical Analysis Latex Fixed volume V Roller the slopes are equal at the contact point R (d-R sin θ, δ+R(1-cos θ)) θ δ d Figure 3.3: Cross-sectional view. The fluid is between two rollers, which have radius R and are located distance d from the center. The fluid follows the circular shape of the roller until a ‘contact point’ at angle θ. (tan θ)2 + 1 d2 h = . (3.10) dx2 R cos θ After selecting reasonable non-dimensional values for parameters involved (such as R and d), equation (3.8) can be solved numerically with MATLAB for various angles θ. See Figure 3.4 for solution and Appendix A for numerical details. Equation 3.8 is a fourth order differential equation with an unknown parameter, the pressure P . Thus, it would require five boundary conditions for a unique numerical solution. To stabilize the system, and prevent the unstable oscillations of the solution, we introduce a second equation to the Z r h2 dx, where y(0) = 0. Then y(r) is the integrated value of system: y = 0 the height squared from the centre to the contact point. The remaining four boundary conditions will require that the height, h(x), and its derivative, h0 (x), are continuous at both contact points. The last remaining unknown parameter is the contact angle, θ in the boundary of the domain. The physics of the problem dictate that the fluid will make the transition at the contact point as smooth as possible, and so the angle θ is chosen such that the second derivatives are continuous as well. There is more than one approach to solving (3.8), and in fact, it can be solved analytically. However, its analytic solution is still presented with a problem of unknown parameters, P and θ. The resulting system of nonlinear equations would still require use of software in order to obtain these 13 3.2. Numerical Analysis 3.5 1.08 0.78 1.3 Roller 3 Height (z−direction) 2.5 2 1.5 1 0.5 0 −3 −2 −1 0 1 Horizontal position (x−direction) 2 3 Figure 3.4: Solutions for different θ’s (as labeled in the legend, in radians). value. We therefore, chose to use a differential equation solver for time efficiency. The obtained shape is used as initial condition in the numerical simulation in the following Section 3.2.3, thus, we are not requiring an analytic form of the solution to the equation (3.8). 3.2.3 Numerical Solution of the PDE The PDE (3.5) can be solved numerically. The details of the numerical computations are given in the Appendix A. The simulation starts with an initial condition, Figure 3.4, satisfying the boundary conditions, specified in the previous section, and computes the shape for the subsequent times. Figure 3.5 shows a sample simulation result with a starting fluid bump held in place by the rollers. As the fluid moves to the sides, the peak flattens out. As that happens, the height begins to decrease more slowly since the fluid motion becomes more gradual and less fluid is leaving from directly beneath the peak. It is worth noting that if one were to consider a constant flux at the center of the fluid spread, then as the fluid propagates the height would remain at the same level, provided the flux is of appropriate magnitude to compensate for the exact amount of fluid that would otherwise 14 3.2. Numerical Analysis 0.9 t=0 t = 20 t = 60 t = 120 t = 220 t = 400 t = 740 t = 2100 t = 5800 t = 30000 0.8 0.7 Height 0.6 0.5 0.4 0.3 0.2 0.1 0 −20 −15 −10 −5 0 5 Horizontal Displacement 10 15 20 Figure 3.5: Numerical simulation of a horizontal fluid bump spread. be left behind. This behaviour resembles a constant flux crack propagation, although the different boundary conditions at the contact line make it a separate problem. Such fractures with constant trailing height have been studied in depth over the past decades [12, 14]. There are several relations of interest in this simulation. First, the propagation of the fluid front. As the fluid spreads out, the fluid film becomes thinner, the height decreases, and the motion slows down. Figure 3.6 shows the behaviour of the front and the height with respect to time on a logarithmic scale, where at larger times this behaviour approaches a linear function. While this cannot be rigorously compared to the experimental results, the overall qualitative characteristics of the behaviour are expected to be the same. When first released the fluid starts spreading rapidly towards a more stable state. Then after a sufficient time the peak is expected to flatten as the fluid spreads out, and the changes in the height and the wavefront are expected to slow down to be linear on a logarithmic scale. This power law behaviour is further examined as a similarity solution in Section 3.4. 15 3.2. Numerical Analysis 0 10 (b) (a) 1.5 Maximum Height of the Fluid Distance traveled by the Fluid 10 Reached Slope −0.197 −1 10 0 10 2 4 10 10 Time Reached Slope 0.211 1.4 10 1.3 10 1.2 10 6 10 0 10 2 4 10 10 6 10 Time Figure 3.6: (a) The height decrease with respect to time. (b) Front propagation with respect to time. 16 3.3. Experimental Results 3.3 Experimental Results The following is a general description of the experiment used to study the behaviour of the fluid on a horizontal plane. For a more technical description of the experimental setup refer to Appendix B. Procedure: Consider fluid covered with an elastic sheet on a horizontal surface. Using two identical rollers with radius R = 1.894 cm, the fluid is collected in the centre (see Figure 3.3). The shape of the fluid does not change in the direction parallel to the rollers, so the cross-sectional view contains all the important information concerning the model. The rollers are removed simultaneously and the spreading of the fluid in both directions is observed. As the fluid is released it immediately starts moving. The thickness at each point decreases as it spreads out, and the motion becomes slower. The points of interest are the shape of the fluid and how fast it moves. In order to compare the experimental results with theory, one must go back to the non-dimensional scaling. As mentioned earlier, one of the scaling parameters has to be determined from the experiment. The parameter of choice is the characteristic height, hc , which is selected to correspond to the initial height of the fluid bump immediately after the release. That is, the characteristic height is the maximum height over the life of the experiment. Figure 3.7 shows a sample frame from which the measurement can be taken. It is the first frame in which the shape of the fluid bump can be observed. The obtained characteristic height value is hc = 0.0059 m and the remaining parameters can be calculated using the experimental physical values given in Appendix B: xc = xc = B ρg 1/4 1/4 0.00533 1400 · 9.81 = 0.02496 m 12µx6c tc = Bh3c (3.11) 12 · 8.5(0.0249)6 0.00533 · (0.0059)3 = 22.5 s. tc = A typical experimental spread is shown in Figure 3.8. The fluid starts out as a relatively sharp peak at the centre and as soon as it is released it rapidly flattens out into a shape that applies less pressure on the elastic 17 3.3. Experimental Results 7 6 Height of the Fluid (mm) 5 4 3 2 1 0 −25 −20 −15 −10 −5 0 5 10 Horizontal Distance (mm) 15 20 25 Figure 3.7: The raw experimental result from the first frame. The shape symmetrized for purposes of numerical analysis. 18 3.3. Experimental Results 9 t=2s t = 23 s t = 43 s t = 103 s t = 203 s 8 7 Height (mm) 6 5 4 3 2 1 0 −50 −40 −30 −20 −10 0 10 20 Horizontal distance (mm) 30 40 50 Figure 3.8: The tallest peak corresponds to the fluid profile 2 seconds after the release and the last profile is three minutes after the release. Three profiles in between are included. sheet. As the boundary moves further away from the centre and the height decreases, the flow gets significantly slower. The first two profiles are 20 seconds apart, while the last two are more than a minute apart. 3.3.1 Experimental Effects of the Bending Stiffness One of the crucial aspects of the experiment is to verify that the behaviour of the fluid is in fact governed by gravity and the bending stiffness of the sheet. It is fairly intuitive that the fluid moves down due to gravity, but the role of the bending stiffness is less pronounced. The main uncertainty is the possibility that the fluid shape is governed by any tension of the elastic sheet accidentally introduced into the setup. Mathematically, the difference between a shape governed by the elastic 19 3.3. Experimental Results plate and surface tension is the difference between solving a sixth order differential equation and fourth order differential equation, respectively. There have been numerous studies, both theoretical and experimental, of the latter and very few of the former. The sixth order problem is more difficult to solve, while the results are qualitatively similar. In both cases the solution yields a spreading fluid and in the case of an inclined plane (as discussed later in Chapter 4) the fluid will have a bump at the moving wavefront. However, the quantitative behaviour can be very different between the two equations. The elastic sheet is resting on top of the fluid, and as the fluid changes shape it has to bend the sheet to match its motion. Meanwhile, the sheet does not transfer any tension to the fluid. A tension in the sheet is only possible if it is pre-stretched when it is fixed on top of the fluid, and if such were the case, the shape of the fluid would be determined by the force from the applied tension. However, unless a very rigorous methodology was used to recreate this imposed tension every time, the experiment would not be reproducible. Figure 3.9 shows how different experimental runs compare. The volume between different experiments does differ in experimental setups but remains of the same order in all of them. The difference can be seen in Figure 3.9, however, the shapes remain similar as the fluid spreads. Two runs of similar starting volume would be indistinguishable from each other, indicating that they can be described by the same model and parameters of that model are statistically the same (unchanged within experimental error). Next, consider the initial shape of the fluid. In this steady state, the shape of the fluid is governed by equation (3.8). Here, the surface tension behaviour would only introduce a second order derivative and a negative sign, so the steady state relation would be, hxx ∝ −h + P. (3.12) Thus, when the fluid is under tension it will behave as a trigonometric function and particularly, since we are considering a symmetric case, a cosine function. A fourth order equation, on the other hand, allows for hyperbolic cosine solutions as well. Figure 3.10 looks at a sample experimental initial shape and compares it with each possible numerical solution. Even before the fluid is allowed to move its initial shape is flatter at the top. While the cosine function approximates the general shape, it fails to 20 3.3. Experimental Results 5 4.5 4 Fluid Height (mm) 3.5 3 2.5 2 1.5 1 0.5 0 −6 −4 −2 0 2 Horizontal Distance (cm) 4 6 Figure 3.9: The fluid shape for several experimental runs. The discrepancy in height is due to different volumes. 21 3.3. Experimental Results 6 5.5 5 Height (mm) 4.5 4 3.5 3 Tension Experimental Bending 2.5 2 1.5 −8 −6 −4 −2 0 2 Horizontal Distance (mm) 4 6 8 Figure 3.10: The initial experimental shape of the fluid (dotted line) approximated by the solutions of the tension equation (dashed) and the bending stiffness equation (solid line). 22 3.3. Experimental Results pick up the slight distortion on the sides, whereas the numerical solution for the bending stiffness problem has a much better agreement with the experimental shape. Another comparison of tension versus bending stiffness is presented in Figure 3.11 where the experimental results are compared to the numerical solution governed by the stiffness of the sheet and to the numerical solution governed by the tension in the sheet. That is, the initial conditions are computed numerically via equation (3.8), Figure 3.4 and simulated as a sixth order PDE model and initial conditions from (3.12) with a fourth order PDE model. The results are compared at a sufficiently long time. The former clearly mimics the shape of the flow more accurately, whereas the latter predicts a much smoother curvature than the experiment. 23 3.3. Experimental Results 0.35 Non−dimensional Height 0.3 0.25 0.2 0.15 0.1 Tension Experimental Bending Stiffness 0.05 0 −10 −8 −6 −4 −2 0 2 Non−dimensional Width 4 6 8 10 Figure 3.11: Comparison of numerical simulations starting with the numerical initial conditions governed by bending stiffness (solid line), and by tension (dashed line) and non-dimensional symmetrized experimental shape (dotted line) after 4 minutes. 24 3.3. Experimental Results 3.3.2 Numerical-Experimental Comparison The shape of the moving fluid is strongly dependent on the initial shape and initial volume of the liquid and in order to be able to compare theory and experiment the starting bump has to be the same in each. One way to achieve this is to take the experimental shape of the fluid and use it as an initial condition for the numerical simulations. Again, the initial shape is taken from the first available frame (Figure 3.7). The experimental and numerical results at various times are compared in Figure 3.12. The starting pre-wetted thickness, δ is selected to correspond to the height of the right-most data point on the experimental curve. Note that such a small value is of the same order as the expected experimental error and thus, it cannot be treated as the true pre-wetted thickness. However, Appendix A shows that changing this thickness by as much as a factor of 10 (from δ = 0.01 to δ = 0.1) has a negligible effect on the results. Observe the last profile has the least correspondence with its numerical simulation. The experimental errors and experimental-numerical discrepancies present at the beginning of the experiment and simulation accumulate as the fluid moves and comparison at later times is not as accurate. The numerical solution of the sixth order problem agrees with the general shape of the fluid. Unfortunately, when the same solution is computed assuming that the tension is governing the motion (solving a fourth order problem), the simulation does not have enough time to evolve to a state where the two solutions are significantly different. That is, the deviation between the two results is less than the experimental error, Figure 3.13. It is not possible to compare the two simulations at a later time because the long time behaviour for the two models is expected to be the same. This is discussed further in Section 3.4. Thus, this comparison neither proves nor disproves the validity of one model over the other. However, as mentioned in the previous section, when the same computation is done using the purely numerical solution, simulated with numerical initial conditions, we can see that the sixth order derivative model mimics the shape of the fluid much more accurately, Figure 3.11. As the fluid spreads, the height of the bump is decreasing and the fluid front is moving forward. This progress can be tracked for both the experimental results and the numerical solution. In addition, this behaviour can be determined through the similarity solution as discussed in Section 3.4. The height of the fluid is simply the maximum point on the curve, which 25 3.3. Experimental Results Time: 1.3 seconds Time: 20 seconds 6 Height (mm) Height (mm) 6 4 2 0 20 40 Width (mm) Time: 72 seconds 60 0 20 40 Width (mm) Time: 208 seconds 60 0 20 40 Width (mm) 60 6 Height (mm) Height (mm) 2 0 0 6 4 2 0 4 0 20 40 Width (mm) 60 4 2 0 Figure 3.12: Comparison of numerical simulation (solid line) and experimental results (dotted line) at different times. Slight discrepancies in the volume of the fluid are due to experimental errors when measuring the fluid height. 26 3.3. Experimental Results 5 Tension Bending Stiffness Experiment 4.5 4 Height (mm) 3.5 3 2.5 2 1.5 1 0.5 0 0 5 10 15 20 25 30 Horizontal Distance (mm) 35 40 45 50 Figure 3.13: Comparison of numerical simulations with bending stiffness (solid line) and with tension (dashed line) at time 9 seconds. Both use experimental steady state shape as initial conditions. 27 Front Position (mm) 3.3. Experimental Results Experimental (Slope: 0.193) Experimental Simulation (Slope: 0.211) 1.5 10 1.4 10 1.3 10 0 1 10 2 10 3 10 10 Time (s) 0.7 Height (mm) 10 0.6 10 Experimental (Slope: −0.189) Experimental Simulation (Slope: −0.192) 0.5 10 −2 10 −1 10 0 1 10 10 2 10 3 10 Time (s) Figure 3.14: Top: Numerical, with experimental initial conditions, simulation (solid line) and experimental (dotted line) front positions with respect to time on a logarithmic scale. Bottom: Numerical, with experimental initial conditions, simulation (solid line) and experimental (dotted line) fluid heights with respect to time on a logarithmic scale. remains near the origin. The fluid front position is the horizontal distance that the fluid has traveled, which is the point at which the height drops close enough to the pre-wetted thickness. Figure 3.14 shows a clear agreement between the numerical simulation, which starts with experimental initial conditions and labeled as ‘experimental simulation’, and the experimental data. Both height and and width change at a similar rate, with a slight difference at later times when the experimental-numerical discrepancies accumulate. 28 3.4. Similarity Solution 3.4 Similarity Solution In the two previous sections we discussed the full numerical solution to the problem and the experimental results. To get a further understanding of the fluid flow, we look for a similarity solution. This is a well known approach to analyzing the long term behaviour of the spread. The analysis presented is very similar to the previous works concerning the spread of the fluid under surface tension [9, 7]. The full governing equation (equation (3.5)) of the flow is ht = h3 hxxxxx + h3 hx x . Far from the fluid front the bending stiffness is less important and the higher order derivatives of h are small. Thus, the time derivative is mostly influenced by the lower order derivative, ht = h3 hx x . Now, assume the similarity solution is of the form x . (3.13) xf (t) here η is the self-similar variable, which measures the relative distance to the wave front, xf (t), as it moves with time. Next, we assume that h0 (t) = Atα , xf (t) = Ctβ and substitute the solution of this form into the governing equation in order to solve for the self-similar exponents, α and β, giving h(x, t) = h0 (t)H(η), η= ∂(h0 H) = h30 H 3 (h0 H)x x , ∂t 2 xf ḣ0 H − h0 H 0 xẋf = 3h40 H 2 (H 0 )2 + h40 H 3 H 00 . (3.14) Here, the 0 indicates the derivative with respect to η and the · indicates the time derivative. The volume is conserved, as is the volume per unit R1 width of the domain. Thus, V = h0 xf 0 Hdη is constant. Differentiating with respect to time, t, gives ḣ0 xf = −h0 ẋf and, consequently, α = −β. Rewriting equation (3.14) to isolate the time dependence, the following is obtained: x2f ḣ0 H xẋf 0 − 3 H = 3H 2 (H 0 )2 + H 3 H 00 . h40 h0 (3.15) 29 3.4. Similarity Solution The explicit time dependence in the terms on the left can be expressed as follows: x2f ḣ0 C 2 α 2β−3α−1 = t , h40 A3 xẋf C 2 β 2β−3α−1 = t η. h30 A3 Since the right hand side of the equation (3.15) has no explicit dependence on time, neither will the left hand side. Hence, either the time exponent will be equal to zero or the left hand side will identically be zero. It can be quickly verified that the latter allows only trivial solutions. Namely, the left hand side requires H = const./η, but the right hand side then evaluates to 5(const.)4 /η 6 , which is non-zero for non-trivial solutions. Therefore, we require the time exponent to be zero, 2β − 3α − 1 = 0, giving the self-similar exponents α = −1/5 and β = 1/5. Then the selfsimilar solution is x . h(x, t) = At−1/5 H Ct1/5 The fluid propagates at t1/5 , and it is thinning at t−1/5 . Once time dependence is removed, the remaining equation is an ODE and can be solved analytically, for − C2 (H + ηH 0 ) = (H 3 H 0 )0 . 5A3 Integrating over the domain of η with condition H 0 (0) = 0 gives ηH(η) + 5A3 H(η)3 H 0 (η) = 0. C2 Solving the differential equation, the η dependence is seen to be r H(η) = 3 3C 2 3C 2 2 C − η , I 5A3 10A3 (3.16) where CI is the constant of integration. Then the complete self-similar solution for h(x, t) is 30 3.4. Similarity Solution r h(x, t) = 3 3C 2 3 x2 CI t−3/5 − . 5 10 t (3.17) This solution relates the height of the fluid to the horizontal distance, x, and the time, t, after its release. We are interested in the maximum height of the bump, hb and the corresponding displacement xf , which is the position of the front for each time t. Ultimately, we would like to compare the behaviour of this solution to the already existing experimental and numerical results from Figure 3.14. In order to do so, consider (3.16), which relates the height to the variable η. From (3.13), we know that at the front position η = 1 and the height should drop to 0, giving r 2 3C 2 3 3C C − , 0= I 5A3 10A3 1 CI = . 2 There is still an unknown front position constant, C, in the solution, which is determined from the volume conservation, where the volume of the fluid, V , is a known parameter. Then volume conservation gives the relationship in terms of the gamma integral r Z 1 p 5 3 3C 3 2 1 − η dη , V = 10 0 #3/5 "r 10 15 Γ(2/3)Γ(5/6) 3 √ C(V ) = V 3/5 . 3 2 3 π 3/2 The front position is approximated by xf = C(V )t1/5 . (3.18) The volume of the fluid in the experiment is known, the constant C can be calculated numerically, giving xf as a function of time. The relation between the similarity wavefront and the experimental/numerical results is shown graphically in Figure 3.15. The similarity solution predicts the long term behaviour to be linear on a logarithmic scale and the figure shows that the curves corresponding to experimental results and both simulations described 31 3.4. Similarity Solution in Sections 3.2.3 and 3.3.2 are approaching the linear behaviour. The final recorded slope of each curve is very close to the predicted magnitude of 0.2. 2 10 Front Position, xf Experimental (Slope: 0.193) Experimental Simulation (Slope: 0.211) Similarity Power Law (Slope: 0.2) Numerical Simulation (Slope: 0.17) 1 10 0 10 1 2 10 10 3 10 Time, t Figure 3.15: Position of the front as a function of time. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. Next, we use (3.17) to calculate how the maximum height at the center of the fluid is decreasing. The centre of the spread corresponds to the coordinate x = 0. That is, the power law between height and time is hb (t) = t−1/5 . Likewise, the relation between decreasing bump height and time is compared in Figure 3.16. 32 3.4. Similarity Solution 0.8 10 0.7 10 0.6 10 0.5 Heigt, hb 10 0.4 10 0.3 10 0.2 10 Experimental (Slope: −0.189) Similarity Power Law (Slope: −0.2) Numerical Simulation (Slope: −0.197) Experimental Simulation (Slope: −0.192) 0.1 10 −2 10 −1 10 0 1 10 10 2 10 3 10 Time, t Figure 3.16: Height of the fluid as a function of time. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. 33 3.4. Similarity Solution 1 Height, hb 10 Experimental (Slope: −1.08) Similarity Power Law (Slope: −1) Numerical Solution (Slope: −1.1) Experimental Simulation (Slope: −1.1) 0 10 1.2 10 1.3 10 1.4 10 Front Position, xf 1.5 10 Figure 3.17: Height of the fluid as a function of the position of fluid the front. Compares the similarity solution (solid straight line) with experimental data (dots), simulation with experimental initial conditions (solid curve) and simulation with numerical initial conditions (dashed curve). The final slope of each curve is stated in the legend. The last comparison is between the height and the corresponding xf value. As the front moves forward, the fluid spreads out, so consider the bump height as a function of the moving front, hb = hb (xf ), and use (3.18) to remove the time dependence from (3.19): Figure 3.17. hb = C(V ) . xf Overall, the slopes correspond to the power law of the similarity solutions. In some cases (like Figure 3.15) the numerical simulation has not fully reached the long time behaviour, but it is clearly approaching the required slope. 34 Chapter 4 Inclined Plane: Numerical Model Chapter 3 discussed the dynamics of the fluid spreading on a horizontal plane. It verified the applicability of the numerical model to the described physical set-up and justified the validity of the experiment. This chapter recalls the more general case of a variable angle, α. The analysis presented is similar in structure, but carries several important qualitative differences. Particularly, we are interested in the instabilities of the fluid wave front, which were not present when the flow was horizontal. 4.1 Non-dimensional Form The general governing equation (2.14) derived in Chapter 2 was h i 1 ∂h 3 4 3 3 = ∇ · Bh ∇∇ h + ρgh cos α∇h − ρgh sin αî . ∂t 12µ (4.1) The non-dimensional scaling now has an extra term (the y-direction), so consider the new co-ordinate transformation to be (x̄, ȳ, h̄, t̄) = (x/xc , y/xc , h/hc , t/tc ). Here, the characteristic in-plane length is taken to be xc , that is, both x and y coordinates have the same scaling. In order to determine a reasonable relationship between the characteristic parameters, one balances the terms of equation (4.1). We first balance the ‘bending stiffness’ term (the first term containing B) with the downwards gravity term (the third term of the right hand side), Bh3 d d4 h ∼ ρgh3 sin α. dx dx4 (4.2) 35 4.1. Non-dimensional Form Since h is of the order hc , we can substitute hc for h in the equation. Also, the derivatives with respect to x becomes 1 d d = , dx xc dx̄ and similarly for the derivative with respect to y. So equation (4.2) yields the appropriate choice for xc : xc = Bhc ρg sin α 1/5 . (4.3) Next, we balance the left hand side of the equation (4.1) with the first term (the ‘bending stiffness’ term) of the right hand side. That is, we require that the time derivative is of the same order as the right hand side, 4 1 d ∂h 3 d d ∼ Bh h . (4.4) ∂t 12µ dx dx dx4 The dimensional scaling for time is then derived by using equation (4.3) to simplify (4.4), giving, tc = 12µxc . ρg sin αh2c (4.5) The two equations (4.3) and (4.5) describe characteristic length and time, respectively, which are the two parameters required to find the velocity. Thus, define the characteristic velocity to be ρg sin αh2c xc = . (4.6) tc 12µ The non-dimensional scalings given in (4.5) and (4.3) are used to reduce equation (4.1) to its non-dimensional form. Rewriting the equation without the bars for simplicity, this gives U= 3 3 ∂h B tc h3c ρg tc h3c 4 = ∇ · h ∇∇ h + cos α∇ · h ∇h ∂t 12µ x6c 12µ x2c ρg tc h2c ∂h3 − sin α 12µ xc ∂x 3 3 ∂h3 4 = Gb ∇ · h ∇∇ h + Gd ∇ · h ∇h − Gi . ∂x (4.7) (4.8) 36 4.2. The Propagating Wave Front Straightforward algebra shows that the coefficients Gb and Gi reduce to unity, which is expected since the scalings were chosen to do so. The only remaining coefficient is Gd , which can be written in various forms, Gd ρg tc h3c hc cot α = hc = cos α = 2 12µ xc xc 1/5 12µ 2 xc = h cot α. B c tc ρg sin α Bhc 1/5 cot α This is the third linearly independent relation between the characteristic parameters. There is no unique way to express the non-dimensional constant Gd = D and the form chosen depends on the known values of the characteristic scalings. It also depends on a number of parameters, so it does not represent the effect of a specific physical quantity, but the overall effect of the physical properties involved. The final non-dimensional form of the governing equation is 3 3 ∂h3 ∂h 4 = ∇ · h ∇∇ h + D∇ · h ∇h − ∂t ∂x 1/5 12µ 2 where D = hU cot α. B c (4.9) This is not the only way to non-dimensionalize the problem. Alternative scaling can be used by balancing different terms of the equation (4.1), in which case a different non-dimensional constant would have a different role in the equation. One such possible scaling is considered later in Section 5.3. 4.2 The Propagating Wave Front Equation (4.9) is a two dimensional partial differential equation, which cannot be solved analytically. Even a numerical solution is nontrivial to obtain. The height of the fluid, h(x, y, t), depends on the x-position, where the xdirection is down the incline, and the y-position, where the y-direction is the direction perpendicular to the flow of the fluid. The downwards motion will evolve more rapidly than the changes to the fluid shape in the y-direction, which is the direction along which the instabilities are formed. We first look at the overall motion of the fluid by using the one-dimensional governing equation to model the wave front, 37 4.2. The Propagating Wave Front ht = h3 hxxxxx + D(α) h3 hx x x − (h3 )x . (4.10) The standard approach in this case is to consider the traveling wave solution to obtain a solution to the steady wave front propagation. Introduce a variable substitution ξ = x − U t and redefine f (ξ) = h(x, t). This turns the governing equation (4.10) into an ODE, U f + f 3 fξξξξξ + Df 3 fξ − f 3 = c, (4.11) where c is the integration constant and U is the speed of the propagating wave. These constants can be determined using the proper boundary conditions. As mentioned in Section 3.2, a pre-wetted layer of thickness δ is introduced to resolve the numerical difficulties at the contact point. Also, since the variables are dimensionless, we choose the non-dimensional height h to be unity behind the wave front. Adding the fact that far from the wavefront the fluid will tend to flatten out, the boundary conditions on a domain of length L are h(0, t) = 1, hx (0, t) = hxx (0, t) = 0, h(L, t) = δ, hx (L, t) = hxx (L, t) = 0. The constant height at the left boundary is equivalent to having a constant fluid source so that the volume is conserved at the propagating front. However, realistically the total volume has to be conserved, so the fluid level behind the moving front cannot be fixed. Rather, its height is decreasing very slowly. Over sufficiently small periods of time the change in height due to volume loss is below experimental error. This makes the approximation valid for our computational purposes, particularly for calculating the stability of the wave front, that is, the stability in the y-direction. The constants c and U can be calculated in terms of δ from the boundary conditions. Particularly, at the boundary conditions the derivative of (4.11) vanish and the height becomes constant, so the equations at the left hand and right hand boundaries reduce to U (1) − (1)3 = c, U (δ) − (δ)3 = c, respectively. The constants c and U are then given by 38 4.2. The Propagating Wave Front U = 1 + δ + δ2, c = δ + δ2. (4.12) (4.13) From the equations it can be seen that for increasing δ the speed of the flow increases as well. As the pre-wetted thickness gets closer to the height at the left hand boundary, which is taken to be unity in this non-dimensional case, the flow moves faster. Also, as δ → 0, U → 1 and becomes independent of δ. The shape of the fluid can be computed from the equation (4.11), and since this equation is a one-dimensional ODE, we can use a built-in MATLAB solver to obtain the solution. The details of this implementation are in Appendix A. The numerical solution for different parameters is given in Figure 4.1. The steady propagating shape captures the bump at the wavefront and the small oscillations around it. These attributes are present in the physical problem even as the fluid flow becomes thinner, but certainly, the height of the bump will decrease as the height of the fluid decreases (Section 4.3.4). As the pre-wetted thickness, δ, increases, the maximum height of the wave decreases and the oscillations become less pronounced. Conversely, for smaller δ’s the peak grows larger and the limit predicts that it would grow infinitely large if δ is taken to zero. This behaviour is similar to the case of surface tension, which is expected since the mathematical analysis of the two problems shares several key features. Thus, the singularity behaviour is similar to that in available literature [15]. An infinitely large peak is not physically reasonable and indicates that the model breaks for a sufficiently small δ. Also, when the pre-wetted thickness increases, it overwhelms the height of the peak, so a sufficiently large δ obscures the behaviour we are interested in. Typical values of pre-wetted thickness are 0.01-0.1 [15, 7], so we use δ = 0.05 for our numerical simulations. Another similarity between the two problems is the effect of the dimensionless constant D on the height of the peak. Higher values of D give a more stable flow. After a certain threshold, when D is greater than some critical value Dc , the peak will disappear altogether and the flow will be stable. [16] The horizontal flow described earlier, α = 0, did not have a peak that forms at the wavefront of the flow, and thus, the flow did not exhibit instabilities in the y-direction. See the analysis in Section 4.3.2. 39 4.2. The Propagating Wave Front MATLAB Solution for D = 1 and Variable δ 2 δ = 0.05 δ = 0.1 δ = 0.2 Height 1.5 1 0.5 0 −5 0 5 Horizontal (x) 10 15 MATLAB Solution for δ = 0.05 and Variable D 2 D=0 D=1 D=4 Height 1.5 1 0.5 0 −5 0 5 Horizontal (x) 10 15 Figure 4.1: Top: Different δ’s. Bottom: Different D values. 40 4.3. Time Dependent Numerical Analysis One final aspect of the propagating wave front to note are the oscillations present in front of the fluid and immediately behind the peak. The elastic sheet is most deformed in this region, where the height changes most rapidly, and as it tries to recover to its neutral flat state it overshoots the equilibrium and then approaches it in an oscillatory manner. 4.3 Time Dependent Numerical Analysis An alternative way to obtain the steady state solution is to solve the PDE (4.10) for longer times, when the solution becomes a steady traveling wave. This method requires considerably more work and involves longer computations. However, it is much more stable with respect to extreme values of δ and D and provides a solution that can change with time. Thus, it can be extended to the cases when the traveling wave shape is no longer fixed, such as a flow with a constant volume. 4.3.1 Solution with Constant Flux Boundary Conditions The numerical methods used to solve the one-dimensional governing equation (4.10) are similar to the ones implemented in Chapter 3. The main differences are outlined in Appendix A. The two major changes, which occur when going from the horizontal flow to the motion on an incline is that the current governing equation, ht = h3 hxxxxx + D(α) h3 hx x x − (h3 )x , (4.14) has an additional gravity term and a dimensionless constant. The numerical scheme also uses Finite Difference for space discretization and an implicit time-stepping method for the time derivative. The numerical scheme starts with an initial condition and evolves with time. The convergence to the propagating front solution is most rapid when the initial condition is a smooth function. Typically, exponential or hyperbolic functions work well. The initial conditions require a function that starts out flat at a constant height, then drops to the constant pre-wetted layer. A candidate for a function satisfying these properties is the hyperbolic tangent, 41 4.3. Time Dependent Numerical Analysis PDE Solution for D = 1 and Variable δ at Time 16 2 δ = 0.01 δ = 0.05 δ = 0.1 Height 1.5 1 0.5 0 0 5 10 15 20 25 Horizontal Distance (x) 30 35 40 PDE Solution for δ = 0.05 and Variable D at Time 16 2 D=0 D=1 D=2 Height 1.5 1 0.5 0 0 5 10 15 20 25 Horizontal Distance (x) 30 35 40 Figure 4.2: Top: Different δ values. Bottom: Different D values. HI = 1 1 − tanh(2x − 20). 2(1 + δ) 2(1 − δ) (4.15) A different function may be used as an initial condition, however, if the chosen function is not sufficiently smooth then the solution may take longer to converge. The numerical scheme may require a smaller time step, because if the time step is too large, the solution from the previous iteration will not suffice as an initial guess, so Newton’s method will not converge. There is no unique way to choose the initial condition, since for the most optimal computation time the starting point will have to be the exact solution. Thus, the hyperbolic tangent is chosen because it is known as a fairly smooth function and the simulation converges in reasonable time. The full numerical solution is given in Figure 4.2. 42 4.3. Time Dependent Numerical Analysis Here, the shape of the fluid varies with values of δ and D in the same way as does the ODE (4.11) solution. This is expected since this is an alternative approach to solving the same problem. One extra observation from Figure 4.2 (Top) is that this flow for larger δ values is faster. This is consistent with equation (4.12), which indicates that the constant speed U will increase slightly with δ. 4.3.2 Stability Analysis Up until now, the one dimensional numerical solution only describes the one-dimensional component profile of the flowing fluid. In order to study the behaviour of the instabilities, the y-dependence has to be brought back into the equation. The vector form of the full non-dimensional governing equation is ∂h3 ∂h = ∇ · h3 ∇∇4 h + D(α)∇ · h3 ∇h − , ∂t ∂x 1/5 12µ 2 hU cot α, where D(α) = B c and ∇ = (∂x , ∂y ). (4.16) Note that the left hand side of the equation (4.16) is a scalar, since it is a derivative with respect to time of a scalar function. Thus, the right hand side can be simplified into its scalar form by expanding the vector products. Take D = D(α) for simplicity. Then ∂h = h3 (hxxxxxx + 3hxxxxyy + 3hxxyyyy + hyyyyyy ) ∂t +(h3 )x [hxxxxx + 2hxxxyy + hxyyyy ] + (h3 )y [hxxxxy + 2hxxyyy + hyyyyy ] +D (h3 )x hx + h3 (hxx + hyy ) + (h3 )y hy − (h3 )x . Now consider small perturbations at the base of the moving wavefront of the form h(ξ, y, t) = h0 (ξ) + h1 (ξ, y, t) where << 1 and h1 ∼ O(1). In this case, h0 (ξ) is the unperturbed solution, which is later shown to be the solution of equation (4.11). These perturbations are expected to evolve in time and to grow larger when the fluid flow is unstable. We will use the change of variables ξ = x − U t to transform the solution into a moving reference frame defined by (ξ, y). Then the time derivative becomes 43 4.3. Time Dependent Numerical Analysis ∂h ∂h0 ∂ξ ∂h1 ∂h1 ∂ξ = + + ∂t ∂ξ ∂t ∂t ∂ξ ∂t = −U h0ξ + (h1t − U h1ξ ) . (4.17) Since h0 (ξ) is independent of y, several of the terms reduce to 0, which yields the final version of the governing equation in a moving frame as follows: −U h0ξ + (h1t − U h1ξ ) = h30 + 3h20 h1 h0ξξξξξξ + h30 [h1ξξξξξξ + 3h1ξξξξyy + 3h1ξξyyyy + h1yyyyyy ] + h30 + 3h20 h1 ξ h0ξξξξξ + (h30 )ξ [h1ξξξξξ + 2h1ξξξyy + h1ξyyyy ] h +D h30 + 3h20 h1 ξ h0ξ + (h30 )ξ h1ξ + h30 + 3h20 h1 h0ξξ + h30 (h1ξξ + h1yy ) − h30 + 3h20 h1 ξ . (4.18) The leading order terms of equation (4.18) are of the order O(0 ) and they give the leading order behaviour as −U h0ξ = h30 h0ξξξξξ + D h30 h0ξ ξ − h30 ξ ξ , (4.19) which can be integrated with respect to ξ to recover equation (4.11), U h0 + h30 h0ξξξξξ + Dh30 h0ξ − h30 = c, (4.20) with U = 1 + δ + δ 2 and c = δ + δ 2 . Thus, there is no new information given at the leading order, but a confirmation of already established behaviour. Next, consider the linear equation for h1 , obtained from the order O() terms: h1t − U h1ξ = 3h20 h1 h0ξξξξξξ + h30 [h1ξξξξξξ + 3h1ξξξξyy + 3h1ξξyyyy + h1yyyyyy ] + 3h20 h1 ξ h0ξξξξξ + (h30 )ξ [h1ξξξξξ + 2h1ξξξyy + h1ξyyyy ] h i 2 3 2 3 +D 3h0 h1 ξ h0ξ + (h0 )ξ h1ξ + 3h0 h1 h0ξξ + h0 (h1ξξ + h1yy ) − 3h20 h1 ξ . (4.21) 44 4.3. Time Dependent Numerical Analysis In order to extract the dependence of the growth rate of the instabilities on the wavelength, consider the perturbations h1 to be in their Fourier Integral form on the infinite domain [−∞, ∞]. Mathematically, this can be expressed in the form Z ∞ g(ξ, q, t)eiqy dq, (4.22) h1 (ξ, y, t) = −∞ where q is the wavenumber of the perturbation, which is related to the wavelength λ by q = 2π/λ, and g(ξ, q, t) is a function containing the spacial and time dependence. Since equation (4.21) is linear in h1 , superposition allows us to consider each q separately. The contribution from each q can be found by substituting h1 (ξ, y, t) with g(ξ, q, t)eiqy rather than the integral (4.22) into equation (4.21), giving, ∂g = ∂t 6 3∂ g h0 6 ∂ξ 4 3 2 3 ∂ g + − − 2q + 3q h0 + Dh0 ∂ξ 2 ∂g U + 3h20 h0ξξξξξ + q 4 (h30 )ξ + 3Dh20 h0ξ + D(h30 )ξ − 3h20 ∂ξ h i + (3h20 h0ξξξξξ )ξ − q 6 h30 + D 3h20 ξ h0ξ + 3Dh20 h0ξξ − q 2 Dh30 − 3h20 ξ g. 5 3 ∂ g (h0 )ξ 5 ∂ξ 4 2 3∂ g 3q h0 4 ∂ξ 2 3 3 ∂ g (h0 )ξ 3 ∂ξ This form of the equation can be simplified using the ODE (4.20) to remove higher order derivatives, ∂g = ∂t ∂ 2g ∂ 6g ∂ 5g ∂ 4g ∂ 3g h30 6 + (h30 )ξ 5 − 3q 2 h30 4 − 2q 2 (h30 )ξ 3 + 3q 4 h30 + Dh30 ∂ξ ∂ξ ∂ξ ∂ξ ∂ξ 2 3c ∂g 3c − 2U + q 4 (h30 )ξ + D(h30 )ξ + −q 6 h30 − 2 h0ξ − q 2 Dh30 g, h0 ∂ξ h0 or written in the operator form, ∂g = Lg. ∂t (4.23) 45 4.3. Time Dependent Numerical Analysis Since the equation is homogeneous, we can assume that time dependence is exponential. Therefore, consider g to be of the form g(ξ, q, t) = φ(ξ)eσ(q)t , σ(q) = growth rate Lφ = σφ (4.24) (4.25) Thus, by solving this eigenvalue problem, the growth rate, σ, can be found as a function of the wavenumber, q, or consequently, as a function of the wavelength of the perturbations. When the value of the growth rate is positive, σ > 0, any perturbation, no matter how small, will grow exponentially and develop patterns, exhibiting the unstable behaviour. The eigenvalue equation (4.25) has a similar structure to equation (4.11), which provides h0 . Therefore, it can be solved numerically using a similar procedure to that given in Appendix A. Figure 4.3 presents numerical solutions with different pre-wetted thicknesses, δ, and with D = 1. 46 4.3. Time Dependent Numerical Analysis Initial Condition (Solution of h0) 2 δ = 0.2 Height 1.5 δ = 0.1 δ = 0.05 1 0.5 0 −5 0 5 x q vs σ 10 15 0.15 σ 0.1 0.05 0 −0.05 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 q Figure 4.3: Top: Leading order solution used as h0 . Bottom: Solution of the eigenvalue σ as a function of q. The mode that will grow preferentially is the q corresponding to the maximum value of σ(q). The instabilities with the fastest growth rate are expected to dominate in the system, so the fluid will have a default preference to form instabilities at the wavelength corresponding to the maximum growth rate. From the figure we can also conclude that as δ increases, the flow becomes more stable because the range of positive σ values shrinks. Hence, there are fewer unstable modes. This is consistent with the decreasing bump height for higher pre-wetted values. The relation between δ and the maximum growth rate is given in Figure 4.4. Again, the figure indicates that the instabilities grow faster as δ decreases. So from the theoretical point of view, the growth rate would be infinitely large as the pre-wetted thickness drops to zero, which is clearly not physically reasonable. 47 4.3. Time Dependent Numerical Analysis 0.16 Maximum Growth Rate 0.14 0.12 0.1 0.08 0.06 0.04 0.02 −3.4 −3.2 −3 −2.8 −2.6 −2.4 log(δ) −2.2 −2 −1.8 −1.6 Figure 4.4: Growth rate dependence on δ value. The limit of small δ causes an increase in peak height and faster growing instabilities. Increasing non-dimensional constant D, however, causes the maximum peak height to decrease. Thus, the flow becomes more stable and the growth rate of the instabilities decreases as well. This is consistent with the surface tension results [16]. Eventually, we expect for sufficiently large D the flow becomes completely stable. That is, growth rate is always negative and the curve lies completely below the x-axis. 4.3.3 Stability Analysis: Small Wavenumber Limit q1 The small wavenumber limit approximates the behaviour of the growth rate at very low q values, that is, very long wavelengths. Consider asymptotic expansions of the growth rate, σ, and the eigenfunction, φ. Since only even powers of q are present in (4.25), all the odd power coefficients reduce to 0, and the expansions for φ and σ should only contain even powers, 48 4.3. Time Dependent Numerical Analysis φ = φ0 + q 2 φ1 + q 4 φ2 + . . . σ = σ0 + q 2 σ1 + q 4 σ2 + . . . (4.26) (4.27) The terms with lower orders of q are expected to contribute the most to the expansion. Take the first two terms of the expansion, and drop all orders higher than q 2 : φ0 σ0 + q 2 (φ1 σ0 + σ1 φ0 ) = 6 ∂5 3 3 ∂ h0 6 + (h0 )ξ 5 φ0 + q 2 φ1 ∂ξ ∂ξ 4 ∂3 2 3 ∂ 2 3 +q −3h0 4 − 2q (h0 )ξ 3 φ0 + q 2 φ1 + Dh30 φ0 + q 2 φ1 ξξ ∂ξ ∂ξ 3c − 2U + D(h30 )ξ φ0 + q 2 φ1 ξ + h0 3c − 2 h0ξ φ0 + q 2 φ1 − q 2 Dh30 φ0 . h0 Then solve for the q 0 leading order term to obtain an equation for φ0 and σ0 : φ0 σ0 5 3c 3 3 ∂ − 2U φ0 . = h0 5 φ0 + D h0 φ0ξ ξ + ∂ξ h0 ξ ξ (4.28) In order to gain insight on the behaviour of the instabilities, we superimpose a small perturbation around the base of the solution in the moving (ξ, y) frame. The position of the starting wavefront line would move from ξ = 0 to ξ = ξb = −eiqy+σt , which is the magnitude of the perturbation. Since ξb is small, linearizing this boundary condition about ξ, we have 0 = h(ξb , y, t) = h0 (ξb ) + h1 (ξb , y, t) + O(2 ) = −h0ξ (0)eıqy+σt + φ0 (0)eıqy+σt + O(2 ). (4.29) In order to satisfy equation (4.29), the order coefficient should be 0. This is satisfied when φ0 (0) = h0ξ (0). Thus, a reasonable approximation of the eigenfunction is φ0 = h0ξ . Substituting this into equation (4.28) we get 49 4.3. Time Dependent Numerical Analysis φ0 σ0 3c 3 3 = h0 h0ξξξξξξ + Dh0 h0ξξ + − 2U h0ξ . h0 ξ (4.30) The next major simplification step is differentiating the first order ODE (4.11) with respect to ξ to show that the right hand side is identically 0, h30 h0ξξξξξξ + Dh30 h0ξξ c + 3 − 2U h0ξ = 0. h0 Hence, σ0 φ0 = [0]ξ . So σ0 = 0. That is, there is no zeroth order approximation to the growth rate and in order to approximate the long wavelength behaviour we have to consider higher order terms. The next order is q 2 , σ1 φ0 = h30 φ1ξξξξξ −Dh30 φ0 + − Dh30 φ1ξ h30 φ0ξξξξ . + 3c − 2U h0 φ1 − 2h30 φ0ξξξ ξ (4.31) Again, this equation is not a trivial one to solve analytically, so we have to take advantage of the known relationships and boundary conditions in order to simplify the result to a more useful form. To solve for σ1 integrate the equation with respect to ξ over the infinite domain: ∞ 3c 3 3 3 − 2U φ1 − 2h0 φ0ξξξ σ1 φ0 dξ = h0 φ1ξξξξξ + Dh0 φ1ξ + h0 −∞ −∞ Z ∞ + −Dh30 φ0 − h30 φ0ξξξξ dξ Z ∞ −∞ As ξ → −∞ and ξ → ∞ the height of the fluid flattens out so the function φ so all of its derivatives are 0. Hence, the terms in the square brackets do not contribute to the integral, i.e. they vanish, leaving Z ∞ Z ∞ σ1 φ0 dξ = −Dh30 φ0 − h30 φ0ξξξξ dξ. (4.32) −∞ −∞ 50 4.3. Time Dependent Numerical Analysis The last few simplifications are obtained by substituting the known approximation φ0 = h0ξ , taking into account that σ1 is independent of ξ and using the ODE (4.11) to replace the higher order derivatives with dependence on δ. We obtain R∞ σ1 = (h30 − U h0 + c) dξ . h0 (−∞) − h0 (∞) −∞ In this derivation we can consider the actual physical problem, since the singularity caused by the contact line does not affect the derivation. Thus, assume that as ξ → −∞, the height will approach unity and as (ξ → ∞), the height will approach zero. We obtain Z ∞ h0 h20 − 1 dξ. σ1 = (4.33) −∞ Thus, the growth rate can be approximated by a parabola, σ = σ1 q 2 , and for small q, the stable flows would be the ones that start out with the parabola opening down (negative growth rate) and unstable ones would have the parabola opening up (positive growth rate). This depends on the sign of σ1 , which can have positive values when h0 > 1 for an interval large enough to contribute to the integral. Hence, the instabilities will form when the fluid forms a large enough bump at the moving front. Figure 4.5 shows how the small limit compares to the growth rate calculation from the eigenvalue problem. To see under which conditions the flow becomes unstable, one considers the small q approximation of the growth rate for various values of D. Solving the zeroth order equation to obtain a solution for h0 at each D value and then using equation (4.33) it is possible to see for which D the growth rate drops completely below the horizontal axis. Figure 4.6 shows that this threshold quantity is D ≈ 3.3. 51 4.3. Time Dependent Numerical Analysis 0.6 Dimensionless Growth Rate, σ 0.4 0.2 0 −0.2 −0.4 −0.6 growth rate at 60o small q approximation −0.8 −1 0 0.1 0.2 0.3 0.4 0.5 0.6 Wavenumber (q) 0.7 0.8 0.9 1 Figure 4.5: A comparison of the small order q approximation with the full solution for growth rate, σ(q), at inclination of α = 60o (D=0.1 and δ = 0.05). 52 4.3. Time Dependent Numerical Analysis −3 2 x 10 D=3 1.5 Growth Rate, σ 1 0.5 D = 3.25 D = 3.3 0 D = 3.35 −0.5 D = 3.5 −1 −1.5 0 0.01 0.02 0.03 0.04 0.05 0.06 Wavenumber, (q) 0.07 0.08 0.09 0.1 Figure 4.6: The approximated growth rate for small q values for different values of D. 53 4.3. Time Dependent Numerical Analysis 4.3.4 Solution with Constant Volume Boundary Conditions Until now, the numerical solution and the stability analysis based on it have considered a model where fluid is being added at a constant rate, which allows for a constant traveling wave solution. This solution is important in order to be able to calculate the growth rate as a function of the unstable wavelength, however, it is not in perfect correspondence with a realistic scenario where the volume of the fluid must be conserved. It is later shown experimentally, in Chapter 5, that close to the moving front this constant flux approximation is reasonable and the stability calculations can still be used for the constant volume case. Nevertheless, it is worthwhile to solve the PDE (4.10) with different initial and boundary conditions to model the fluid flow with conserved volume. Due to constant volume loss at the wavefront, there is no steady state solution which can be used for stability analysis, so all the information we have about this case is that which can be obtained from the numerical solution of the PDE. The initial shape of the fluid is a combination of several exponential functions, since they tend to be smooth enough for the method to easily converge, HI = x2 1 1 − tanh(2x − 20) − (1 − δ)e− 10 . 2(1 + δ) 2(1 − δ) (4.34) A different smooth function may be used as the initial condition as well, for example, a more symmetric curve. However, in that case the starting shape takes longer to converge to the right solution, which significantly increases the overall run time. The changes in boundary conditions are fairly straightforward. Instead of maintaining a constant height of unity on the left hand side boundary, the volume conservation requires a constant height of δ, to account for the pre-wetted layer, and vanishing derivatives. Thus, the boundary conditions become symmetric. The full solution is given in Figure 4.7. As expected, the thickness of the fluid decreases as more fluid gets left behind; the speed of the wavefront decreases as well and at sufficiently long times the remaining fluid stretches into a triangular tail. Yet, the oscillations behind and ahead of the moving front are still present and immediately behind the capillary ridge the fluid oscillates at a constant level before it starts to steadily decrease. 54 4.3. Time Dependent Numerical Analysis Constant Volume: D = 1, δ = 0.05 1.4 1.2 t=6 t=0 1 Height t = 18 0.8 0.6 0.4 t = 42 0.2 0 0 5 10 15 20 25 Horizontal Distance 30 35 40 Figure 4.7: Numerical solution of the PDE (4.10) with the constant volume boundary conditions. 55 4.3. Time Dependent Numerical Analysis A number of variables can be tracked while the solution is changing, such as the peak height, how fast the position of the peak is moving and the slope of the trailing fluid. These values remained constant for the solution assuming a constant flux and they do not change significantly over short periods of time in the constant volume case. Figure 4.8 shows a sample simulation at an angle of 60o . 1.4 1.2 Height 1 0.8 0.6 0.4 0.2 0 0 5 10 15 20 25 Horizontal Distance 30 35 40 Figure 4.8: Tracking the shape of the wave front over a long time (up to t=50). This tracks the height and the position of the peak with time. The behaviour of these parameters can be converted to their dimensional form via the dimensional scaling equations (4.3) and (4.5) and compared with the experimental results in the next Chapter. Plotting the position of the peak versus time on a log-log scale provides insight as to how the fluid is moving (see Figure 4.9). The front behaviour becomes linear with time on a logarithmic scale, indicating that it behaves as a power law, propagating close to O(t1/3 ). 56 4.3. Time Dependent Numerical Analysis 6 5.8 Log(Peak Position) 5.6 5.4 5.2 5 4.8 Numerical Points Linear Fit 0.3826t+4.4379 4.6 4.4 0 0.5 1 1.5 2 Log(time) 2.5 3 3.5 4 Figure 4.9: Position of the wave front against time (dimensionless form). 57 Chapter 5 The Elastic Experiment The complete description of the experimental procedure is described in Appendix B. The general structure of the experiment is similar to the experimental setup of the flow on a horizontal plane. However, the fluid now starts on one side of the base, which is then inclined at a specific angle α. The fluid forms instabilities as it flows down the plane. Their shape, specifically growth rate and wavelength, is an additional aspect of the motion that can be compared with the theory. Overall, the technical aspects of the comparison are similar to those in Chapter 3, with several new parameters. 5.1 Dimensional Variables The numerical equations are reduced to their non-dimensional form prior to solving, so the obtained solutions are in terms of the non-dimensional variables. In order to compare theory with experiment, one again needs to scale accordingly. The non-dimensional variables introduced are: 1/5 Bhc , ρg sin α 12µxc , ρg sin αh2c 1/5 12µ 2 xc h cot α, B c tc xc . tc xc = tc = D(α) = U = (5.1) (5.2) (5.3) (5.4) This defines five variables of interest, that is: xc , hc , tc , U and D(α). However, there are only four equations above, so one of the variables has to be determined experimentally. The chosen parameter is the same as before, that is, the characteristic height hc , which now corresponds to the height of the fluid behind the wavefront. The quantitative value for hc is 58 5.1. Dimensional Variables Table 5.1: The values for hc = 3mm Angle xc (cm) D tc (s) U (mm/s) 30 1.877 0.277 9.476 1.981 45 1.751 0.171 6.252 2.801 60 1.682 0.103 4.902 3.431 78 90 1.641 1.634 0.039 0 4.235 4.125 3.875 3.962 determined several seconds after the fluid’s release, so the fluid has time to form the wavefront. The measurement is taken via the calibration procedure described in Appendix B and its typical values range between 2 − 3 mm. Flow patterns at different angles may vary, but the amount of fluid present remains the same and the starting position is similar, thus the same value of hc can be used. Once the value for hc is obtained, the remaining dimensional parameters can be calculated (including the non-dimensional variable D) since all the other variables are constant properties of the materials used. The next section shows how the characteristic height is chosen for a set of experimental runs and Table 5.1 summarizes the calculated constants for that specific value of hc . There are two numerical solutions that were determined in the previous Chapter. The first approach considered a constant flux, which meant that the height behind the wavefront is constant and equal to unity. The second approach considered constant volume, wherein the height of the front and the flat tail section immediately behind it is constantly decreasing. When the fluid level behind the wavefront is always constant the scaling is straightforward. On the other hand, for constant volume the initial shape does not indicate the height of the fluid behind the front and once the fluid starts moving, that height is never constant. In this case, we can determine the proper scaling by running the numerical simulation for several time units to a time when the peak first forms. Figures 4.8 and 4.7 from Section 4.3.4 show that this happens fairly early in the simulation. Over such short time intervals, the decrease in height is negligible and the initial level can be taken as the characteristic height behind the front. Figure 5.1 indicates where the characteristic height is taken on the experimental and theoretical wavefronts. The horizontal distance is scaled by the characteristic length, xc , to be in 59 5.1. Dimensional Variables hc = 0.5853 Numerical Height 1.5 1 0.5 Experimental Height (mm) 0 0 50 100 150 200 250 Horizontal Displacement 300 350 400 hc = 2.987 mm 8 6 4 2 0 0 20 40 60 80 100 Horizontal Displacement (mm) 120 140 Figure 5.1: Top: Numerical solution (δ = 0.05) with the indicated height used as hc = 0.5853. Calculated for time = 2 which is 9.8 seconds when converted to dimensional units for this image, but this conversion depends on tc and may vary for different inclinations. Bottom: Experimental behaviour with corresponding hc = 2.987 at 10 seconds after the release. This time may vary for individual experiments from 10 to 15 seconds . 60 5.2. Contour Data Processing 8 Height (mm) experimental numerical (δ = 0.05) 6 4 2 0 0 20 40 60 80 Horizontal Distance (mm) 100 120 140 5 experimental numerical (δ = 0.05) Height (mm) 4 3 2 1 0 0 20 40 60 80 Horizontal Distance(mm) 100 120 Figure 5.2: The comparison of experimental and numerical results with angle of inclination of 60o . With tc = 4.9s and hc = 3mm. Top: Comparison 10 seconds after the release with constant volume numerical solution. Fluid is gradually transforming into the shape of the traveling wave front. Bottom: Comparison 1 minute after the release with the constant flux numerical solution. Note, the origin of the horizontal scale has been shifted. the same units as the experimental results. The height comparison between the experimental and numerical predictions is given in Figure 5.2. There is definite agreement between the experimental results and those predicted numerically. One can even observe the oscillations at the base of the front and immediately behind it. Moreover, after a sufficiently long time, the constant flux numerical solution can fairly accurately approximate the experimentally measured shape. 5.2 Contour Data Processing As fluid is moving along the plane, the parameters of interest are the growth rate of the instabilities, the wavelength of the instabilities, and the speed 61 5.2. Contour Data Processing of the wavefront. All of these values are determined by looking at the wave front contour line. A typical set of contours is shown in Figure 5.3. The fluid starts out at the bottom as a flat front and as it moves the fingers start to form. Here, peaks are the tips of the fingers and valleys are the space between the fingers. Typically, four to five peaks form each experimental run. 300 130 133 250 120 Vertical Direction (pixels) 100 90 110 80 200 150 70 60 50 40 30 20 100 10 50 0 time = 0 0 100 200 300 400 Horizontal Direction (pixels) 500 600 Figure 5.3: The contours of fluid fronts at inclination of 78o from time 0 to 133 seconds. 5.2.1 Growth Rate The growth rate is a measure of how fast the instability is growing. That is, it measures the rate of increase in length of the actual finger with respect to time. Intuitively, we expect the growth rate to be larger for a steeper angle. When the incline is sharper, gravity exerts stronger force on the fluid and one would expect the unstable front to grow more rapidly. The instability analysis assumed that the initial growth rate had an exponential form. Thus, it can be computed by plotting the length of the forming finger versus time on a logarithmic scale. In order to determine this slope, we go back to the contour plots as shown 62 5.2. Contour Data Processing (a) 45o 2 10 Log(Length) Log(Length) (b) 60o 2 10 0 10 0 10 −2 10 −2 10 0 50 100 Time (s) 150 0 (c) 60o 2 150 10 Log(Length) Log(Length) 100 Time (s) (d) 90o 2 10 0 10 −2 10 50 0 0 10 −2 50 100 Time (s) 150 10 0 50 100 Time (s) 150 Figure 5.4: Measuring growth rates in Figure 5.3. The length of the finger is the difference between the peak position and the valley position with respect to the vertical axis. The initial growth rate is the growth rate sampled at early times when the instability is first starting to form. Figure 5.4 shows examples of how the slope is measured. Ideally, we look for a curve that starts out as a straight line, indicating the initial exponential growth. The instabilities may not start growing immediately after the release of fluid as can be seen in subplots 5.4(a) and 5.4(c). It is also evident that there is some noise at the beginning of the curve. On occasion, due to experimental error, there is too much noise present in the measurements at early times to determine the slope of the straight line and not every measurement can be easily interpreted. In cases when the data is too vague, that graph is omitted. The measured growth rates are summarize in Table 5.2. The units of growth rate are (s−1 ). Hence, one would expect that longer characteristic times would correspond to smaller growth rates (that is, smaller inclination angles). The characteristic times are shown in Table 5.3. The theoretical prediction for growth rates was done in Chapter 4. There is a fundamental difference between the experimental and predicted growth 63 5.2. Contour Data Processing Table 5.2: Growth rates Angle Measured Growth Rates Average Error 30 0.0287 0.0168 0.0265 0.0296 0.0237 0.0267 0.0339 0.0266 0.0053 45 60 78 90 0.0463 0.0392 0.0433 0.0310 0.0417 0.0778 0.0615 0.0711 0.0684 0.0510 0.1015 0.1217 0.0833 0.0931 0.1033 0.1108 0.1305 0.1272 0.0403 0.0697 0.0901 0.1180 0.0058 0.0067 0.0260 0.0130 Table 5.3: Characteristic time Angle Time tc (s) Growth Rate 30 45 60 78 90 9.476 6.252 4.902 4.235 4.125 0.0266 0.0403 0.0697 0.0901 0.1180 64 5.2. Contour Data Processing rates. Namely, the theoretical calculation relies on a steady propagating wave solution, whereas experimentally the height of the wavefront is continuously decreasing. However, from Figure 5.2(b), we can see that for large times the constant flux numerical solution shows closer agreement with the experimental than the numerical solution with constant volume. Hence, although this difference is one of the sources of error, it is still reasonable to expect consistency between the calculated and measured growth rates. The results are compared in Figure 5.5. The theoretical results are found by taking the maximum growth rates and converting them to their dimensional form using the characteristic values from Table 5.1. 0.14 Experimental Numerical Growth Rate (s−1) 0.12 0.1 0.08 0.06 0.04 0.02 20 30 40 50 60 70 Angle (degrees) 80 90 100 Figure 5.5: Growth rates 5.2.2 Wavelength The wavelength is a measurement of how many fingers can form across a given width. That is, it gives the average distance from one fingertip to another. This measurement is very straightforward, since after the fingers form (see Figure 5.3) the position of the peak does not move, only the length of the finger increases with time. To calculate the wavelength numerically, consider the relationship be65 5.2. Contour Data Processing Table 5.4: Characteristic length with numerical wavenumber and wavelength Angle (Degrees) Length xc (cm) Wavenumber q Num. Wavelength (cm) Exp. Wavelength (cm) 30 45 60 78 90 1.88 1.75 1.68 1.64 1.63 0.445 0.450 0.4590 0.467 0.473 26.5 24.4 23.0 22.1 21.6 16.3 16.7 15.5 14.9 14.9 tween the wavenumber and the wavelength, λ = 2π/q, where q is the wavenumber that corresponds to the maximum growth rate for a specific value of D. Since the wavelength is measured along the y-direction, the scaling used to convert it to its dimensional value is xc . The results along with the experimental measurements are summarized in Table 5.4. According to the theoretical prediction, the wavelength is expected to be longer for smaller angles. However, this variation is very gradual and values across a large range of angles are still very close. Experimental results do tend to decrease slightly for higher angles, but unfortunately this difference is of the same order of magnitude as the experimental error. Thus, it is difficult to obtain results which can state with certainty that the wavelength measurements decrease as the inclination angle increases. This result can also be predicted from the fact that the characteristic lengths, xc , show little variation over different values of the angle α. Since the scale used to convert the numerical values to their dimensional form remain of the same order, it is natural to expect the dimensional value to behave similarly. The experimental results confirm this. Hence, the inclination angle α does not play an important role in determining the wavelength. Other parameters, such as viscosity or elastic sheet thickness, that are currently held fixed in the experiment, may have a much greater effect on the wavelength. 5.2.3 Velocity One more measurement, which characterizes the motion, is the speed of the flow. Numerically, the value can be calculated by simulating fluid flow with constant volume, as described in 4.3.4. The speed of the flow is determined by tracking the position of the highest point (the peak of the wave front bump). 66 5.2. Contour Data Processing The speed of the fluid is determined from the distance it travels down the plane in a fixed time, which is measured by recording the position of the fluid along the y-axis at a fixed x-coordinate in the the contour plot (Figure 5.3) at every frame. Figure 5.6 shows a typical plot of this traveled distance and velocity versus time for the moving fluid. The velocity is the derivative of the distance function, however, since the recorded distance is not a smooth curve due to experimental errors, we cannot determine velocity from taking the slope of raw experimental data. Instead, the distance is fitted with a polynomial and velocity can be easily computed. As expected the fluid starts out faster and then its velocity decreases as more and more fluid is left behind the front. Also, the figure shows the velocity profiles for the fluid near the tip of the finger, which clearly travels faster since the finger is also a growing instability, and for the fluid in between the growing instabilities (the valley between the two adjacent fingers). There is a range of velocities at which the fluid can move and for the purposes of comparison, we consider the average velocity, that is, the velocity of the centre of the forming fingers. 67 5.3. Alternative Dimensional Scaling Distance Traveled by the Fluid Distance (pixels) 250 200 150 Finger peak Valley 100 50 0 0 20 40 60 80 100 120 Time (s) Speed from the Polynomial Approximation 140 160 Speed (pixels/s) 3 Speed of the Finger Speed of the Valley 2 1 0 0 20 40 60 80 100 Time (s) 120 140 160 Figure 5.6: Top: The distance traveled by the wave front from some initial level (the initial level can be chosen arbitrarily, since we are interested in the slope). The black solid line indicates the shape of the fitted polynomial. Bottom: The speed of the wave front, as calculated from the derivative of the distance traveled. 5.3 Alternative Dimensional Scaling The scales used to non-dimensionalize the problem are not unique. This section considers an alternative approach, which will balance the left hand side with the second and the last terms on the right hand side of the equation h i ∂h 1 = ∇ · Bh3 ∇∇4 h + ρgh3 cos α∇h − ρgh3 sin αî . ∂t 12µ (5.5) Using these relations we still get an equation for the each of characteristic length, xc , and time, tc , in terms of the characteristic height, hc : 12µ cos α , ρghc sin2 α = hc cot α. tc = (5.6) xc (5.7) 68 5.3. Alternative Dimensional Scaling Table 5.5: The values for hc = 3mm Angle xc (cm) Gb tc (s) U (mm/s) 30 45 0.520 0.300 615 6780 2.623 1.07 1.981 2.801 The main difference is that the dimensionless constant is now in front of the ‘bending stiffness’ term, so the non-dimensional form of the equation is 3 ∂h3 ∂h 3 4 = ∇ · Gb h ∇∇ h + ∇ · h ∇h − , ∂t ∂x B sin4 α . where Gb = ρgh4c cos5 α (5.8) The coefficient in front of the ‘bending stiffness’ term is expected to be large in order for the bending stiffness to be important to the motion of the fluid. Thus, to verify its importance, one considers the magnitude of Gb based on the experimental parameters. The characteristic constants are given in Table 5.5. Notice that the characteristic velocity remains the same as for the original scaling (see Table 5.1), which is expected since the governing equation has not changed. The new non-dimensional constant Gb , however, is now several orders of magnitudes larger than other involved parameters, indicating that the bending stiffness does in fact play an important role in the modeling of the problem. The PDE (5.8) may also be solved numerically to give the same solution as obtained in Chapter 4. While the solutions are identical, this approach foresees several complications. First of all, by balancing the last two terms of 5.8, we are effectively setting the derivative of h to be of the same order as tan α. This implies that this scaling is only effective at low angles and when α increases by a substantial amount Gb becomes unreasonably large. Secondly, working with such large numbers may cause numerical instability. Thus, a more practical scaling is selected for full model analysis. 69 5.4. Discussion of Experimental Errors 5.4 Discussion of Experimental Errors The experimental results are in good agreement with the theoretical model. However, there are several discrepancies due to experimental errors and some simplifications of the model. The most significant experimental errors come from the procedure used to measure the thickness of the fluid. While measuring light intensity reproduces the general shape of the fluid fairly accurately (see Figure 5.2), it still easily affected by any discolouration of the image, whether it is from the equipment or due to the quality of the image. Moreover, these imperfections are present in both the experimental fluid and the wedge used for calibrating the height. Although the measurements are averaged in order to reduce these effects, it is not possible to remove them completely. Thus, the inaccuracies in calibration are only amplified when the actual fluid height is determined, resulting in a non-smooth curve (see Figure 5.1 (bottom)). Another aspect which should be considered is scale of the experiment. Ideally, we would like a large number of finger instabilities to work with, however, with a wavelength on the scale of several decimeters it becomes extremely difficult to work with more than four or five fingers forming at a time. Thus, the fingers furthest from the centre can potentially be affected by the boundaries of the experimental setup. One more source of errors comes from the known parameters. The values of Young’s modulus, latex sheet thickness and fluid viscosity are all subjected to experimental errors as well and consequently have an effect on the calculations of characteristic parameters. There are also factors which do not have a significant effect on the experimental results, such as homogeneity of the lighting (test measurements were taken to ensure that MATLAB Image Editor does not pick up visible dependence on the lamp position) and calibration fluid (the exact same batch of dyed fluid used for the experiment was used for the calibration). Overall, the results show a remarkable agreement between the numerical and experimental shape of the fluid, particularly at early stages. Also, experimental evidence supports the approximation that the fluid height immediately behind the wavefront is close to constant, justifying the numerical approach to the stability analysis. Moreover, the theoretical growth rate accurately predicts the growth rate at lower angles and estimates the correct order of magnitude for higher angles. It is likely that at higher angles the fluid moves much faster, allowing more room for experimental error. 70 5.4. Discussion of Experimental Errors In general, the growth rate is of order 0.01-0.1s−1 . This is significantly slower than the typical behaviour of the same fluid under surface tension, when the instability growth rate ranges from 0.15s−1 to 0.3s−1 [7]. It is expected that an elastic sheet covering the fluid would make it more difficult for the flow to progress, thus when all the other parameters are kept the same, the growth rate would be much slower. The experiment also confirms the theory that higher bump height corresponds to a more unstable flow, that is, a faster growth rate, see Figure 4.3. Thus, decreasing the pre-wetted thickness, δ, or the constant D would make the flow more unstable. The one parameter the theory does not predict as accurately is the wavelength of the instabilities. While the results are still within 30% of each other, theory greatly overestimates the wavelength in comparison to the experimental results. It is plausible that from the time when instabilities are starting to form (the stage at which the growth rate is measured) to the time when they are starting to grow, the fluid loses a significant portion of its volume behind the wavefront, thus, effectively reducing the characteristic height. This sudden volume decrease may cause the fluid to adjust its original wavelength to a more favourable value. Note that equation (5.1) verifies that decreasing the characteristic height, hc , will decrease the characteristic length, xc , which is in turn is proportional to the wavelength of the instabilities. While minor changes do not have a notable effect due to the fifth root, a significant enough change can decrease the theoretical results up to 20%. More experimental work is needed to justify this hypothesis. 71 Chapter 6 Similarity Approach Section 3.4 introduced the method of similarity solution, which examines the asymptotic behaviour of the fluid at long times. The analysis of the full governing equation is similar, however, non-zero inclination introduces an additional gravity term which will dominate at long times. There is a known solution which has been derived for the surface tension problem in a more general case [9]. Since the method is based on taking the leading term of the governing equation, higher order derivatives become less important and the full governing equation (4.10) can be reduced to ∂h + h3 x = 0. (6.1) ∂t This equation is identical to the one describing the fluid moving under tension and the analysis is consistent with literature [9, 7]. Assume a self-similar solution of the form h(x, t) = h0 (t)H(η), where 0 ≤ η = x/xf (t) ≤ 1 measures the distance to the front xf (t). Once the time dependence is isolated, the shape of the fluid film can be considered as a function of η. The time dependence itself can be expressed as a power law with self similar exponents α and β h0 (t) = Atα , xf (t) = Ctβ . (6.2) Now (6.1) can be simplified to h2 h˙0 xf H − ηH 0 + 3 0 H 2 H 0 = 0. h0 ẋf ẋf (6.3) The conservation of volume still holds, and requires Z 1 0 = ḣ0 xf + h0 ẋf H(η)dη. 0 72 Chapter 6. Similarity Approach This allows for further simplifications to (6.3) and gives one of the equations needed to solve for the self similar exponents: Aαtα−1 Ctβ = −Atα Cβtβ−1 , α = −β or β = −α. Note that this condition is the same as in the case of the horizontal flow. The difference between the two cases comes from eliminating the explicit time dependence, A2 t2α h20 = , ẋf Cβtβ−1 1 3 Hence, we expect xf ∝ t1/3 , so the front propagates with the order t1/3 . Similarly, h0 ∝ t−1/3 , so the height of the wave front bump decreases at the rate of t−1/3 . This is consistent with the earlier derivation from Section 4.3.4, shown in Figure 4.9, that the long time fluid front behavior approaches a power law with respect to time. The exponents for an inclined plane are of higher order than the ones derived for horizontal flow. When the fluid is flowing down an incline then it is expected that the gravity will play a more crucial role and the transformation of the fluid shape will happen faster. 2α = β − 1 ⇒ α=− Equation (6.3) can be solved analytically at this point. Consider −H − ηH 0 + 9 A2 2 0 H H = 0. C (6.4) Integrate with respect to η, to get A2 3 H (η) + CI = 0. C with CI being the constant of integration. As before, the height of the fluid at the front position should drop to zero, which implies that CI = 0. This simplifies the equation and where we can solve for H(η) explicitly, r C H(η) = η. (6.5) 3A3 ηH(η) − 3 73 Chapter 6. Similarity Approach Equation (6.5) only gives the dependence of the spatial part of the solution on η. To get the complete solution substitute (6.5) into h(x, t): r x h(x, t) = h0 (t)H(η) = . 3t Again, the solution depends on two variables, so we consider the conservation of volume to reduce the number of unknowns. Specifically, Z 1r C V = AC ηdη, 3A3 0 1/3 27 2 V . (6.6) C= 4 Once this constant is known in terms of volume, the position of the wavefront xf (t) can be found as a function of time, xf (t) = 1/3 27 2 V t 4 (6.7) Figure 6.1 compares all the obtained results, with the experimental front position being the average distance the fluid has traveled, the numerical simulation obtained by solving the PDE as explained in Chapter 4 and the similarity solution given by 6.7. Initially, the experimental behaviour is consistent with the numerical simulation although it tends to move slightly faster than its numerical counterpart. The differences between experiment and theory are expected. Namely, one can refer to the bottom graph of Figure 5.2 in Chapter 5 to see that at later times, the experiment tends to mimic the behaviour of a constant flux, and thus loses volume slower than the theory, resulting in faster motion of the front. The similarity solution, on the other hand, considers the behaviour of the flow at longer times, so the initial behaviour is expected to be drastically different. At later times however, the behaviour is in fact similar, as all the curves approach a line of constant slope on a logarithmic scale. Similarly, we can find the relationship between the height at the wavefront, hf (t), and time, and consequently, the relationship between hf and xf : 74 Chapter 6. Similarity Approach Distance Travelled by the Fluid 2 Front Position 10 Similarity Numerical Simulation Experimental −1 10 0 1 10 10 2 10 Time Figure 6.1: Position of the front with respect to time. Experimental results (dashed line), numerical simulation (solid line) and similarity solution from equation 6.7 (dotted line). 75 Chapter 6. Similarity Approach hf (t) = V 2t hf (xf ) = 1/3 3V 2xf (6.8) (6.9) The height of the fluid is not tracked experimentally with respect to time, and thus Figures 6.2 and 6.3 only compare the numerical simulation with the similarity solution. Figure 6.2 shows that the height approaches the asymptotic behaviour very slowly, but its dependence on wavefront position rapidly becomes linear, Figure 6.3. While the similarity solution predicts that long term dependence among the measured values is linear on a logarithmic scale, it has several shortcomings. The prediction is identical for flows governed by surface tension and those covered by an elastic sheet. Also, the predicted slope of the behaviour is independent of any problem parameters, such as angle of inclination, viscosity or the bending stiffness. Figures 6.1-6.3 indicate the numerical simulation will eventually reach the linear state, however, in each case this occurs at a different time scale and changing the parameters may affect how fast the fluid behaviour approaches the power law behaviour. The approach does not provide insight on when the similarity solutions becomes a valid solution and at what point in time the behaviour can be approximated with a power law. 76 Chapter 6. Similarity Approach 2 10 Numerical Similation Similarity 1 Front Height 10 0 10 −1 10 −2 10 −1 10 0 1 10 10 2 10 3 10 Time Figure 6.2: Fluid height at the front with respect to time. The similarity solution is given by (6.8). 77 Chapter 6. Similarity Approach Numerical Simulation Similarity 0 Front Height 10 −0.1 10 −0.2 10 −0.3 10 0 10 1 10 Front Position 2 10 Figure 6.3: Fluid height at the front with respect to front position. The similarity solution is given by (6.9). 78 Chapter 7 Conclusion 7.1 Summary In this study, we looked at the theoretical modeling and experimental behaviour of thin fluid film flow instabilities between an elastic sheet and a planar surface. We primarily focused on how the behaviour of the fluid is affected by changing the angle of inclination of the plane from horizontal (angle of 0o ) to vertical flow (angle of 90o ). We established that the introduction of an elastic sheet will affect the initial shape of the flow and its progression. Figure 3.10 demonstrates, through comparison of experimental data with numerical simulations, that the initial shape of the fluid governed by bending stiffness is qualitatively different from the shape of a fluid moving under tension. Figure 3.11 verifies that as the fluid progresses with time, its shape is more accurately represented by the bending stiffness model. The qualitative overlap of analysis of the tension model and the bending stiffness model happens when we look at the long term behaviour of the fluid. At long times, the lower order derivatives dominate, so fourth order (tension) and sixth order (bending stiffness) derivative terms become less important, resulting in the identical asymptotic behaviour of the two models. In our analysis we use a similarity solution to predict the power law, which will approximate the solution at long times. We determined that on a horizontal surface the height of the fluid decreases at the rate of t−1/5 and the wavefront of fluid moves forward at the rate of t1/5 , where t is time. The results were in good agreement with numerical simulations as well as the experimental data. Overall, the flow on a horizontal plane did not exhibit any instability on the scale where the analysis and experiment were performed, however, it showed that the experimental results were consistent with the model predictions, both analytical and numerical, and thus, validated the mathematical model introduced in Chapter 2. 79 7.2. Future Work Introducing the inclination angle results in an additional gravity term acting tangent to the direction of the flow, and the nondimensional constant D, into the governing equation. This causes finger-like instabilities to form and we use the mathematical model to predict the growth rate and the wavelength of these instabilities. The predicted growth rate is consistent with the experimental data and is remarkably accurate at lower angles. The predicted wavelength was of the same order of magnitude as the experimental results, to within 30%. It was also found that increasing the angle has the greatest effect on the growth rate, while the wavelength does not change significantly for steeper angles. Nevertheless, faster growth of instabilities for steeper inclination suggests that the flow is more unstable, which is confirmed both theoretically and experimentally. The long term behaviour of the fluid flow on an inclined plane is also investigated using a similarity solution. However, the dominant term in this case is the tangential gravity component, which is zero in the horizontal flow, and the asymptotic behaviour follows a different power law. The height of the fluid is now decreasing on the order of t−1/3 and the wavefront growing at a rate of t1/3 . The inclined wavefronts evolve more rapidly, Ot1/3 , than the horizontal wavefronts, Ot1/5 due to the enhanced gravitational forcing in the direction of the flow. The comparison between analytical, numerical and experimental results indicates that the fluid behaviour is consistent with the subdominance of the higher order derivatives at large times. The work done in this thesis establishes that the introduced mathematical model is applicable to fluid flows with stiff boundary conditions, through the remarkable consistency between theoretical and experimental results. 7.2 Future Work This work mainly focused on how the angle of inclination affected the nondimensional constant D, and consequently, how it affected the front instabilities. It was noted that the growth rate increases with increasing angle, whereas the wavelength did not have a notable dependence on the angle. There are other physical parameters in the problem which can be varied and the measured dependence compared with theory. Examples include the properties of the fluid, such as viscosity µ or density ρ, and the properties of the elastic sheet, such as the thickness of the sheet. Unfortunately, the 80 7.2. Future Work current experiment cannot accommodate many of these variations, such as varying the thickness of the elastic sheet, since the measured results are not sufficiently precise to see the correlation and introducing a thicker film may introduce additional experimental errors when light intensity is used to measure the height. The relationship between the wavelength of instabilities and the growth rate indicates that the maximum growth rate, with its corresponding wavelength, will dominate and thus be preferred by the fluid front. However, an interesting experiment would be to examine how the behaviour of the fluid changes when a specific wavelength is prescribed. That is, would the fluid front instabilities continue to grow with that wavelength at a slower rate, or would they attempt to readjust the the wavelength to match the preferred scenario. Current experimental setup cannot resolve this question. Due to the experimental error, the measured growth rates would not be statistically different from each other. Although, one way to implement this experimentally would be to enforce small amplitude perturbations to the fluid wavefront line prior to the release of the fluid. Currently, the wavefront is kept completely flat and instabilities form at the favourable wavelength. The study of instability formation can also be extended to incorporate other physical phenomena and their effect on the flow. For example, thin films are known for their wrinkling morphology [6, 13], and while the wrinkling was not a significant factor in this work, a modified study, perhaps on a smaller scale, would be of interest. The only wrinkling present in the current experiment occurred at the very late stages of fluid progression, after the instability had fully formed and all the data had been collected. The presented model is a stepping stone between the fluid flow with a free surface and the fluid flow completely enclosed in a solid. The presence of the elastic sheet constrains the flow, yet, it does not go the extreme where the fluid has to break its own path in order to progress. Thus, while there are a number of experimental studies which can expand on the comparison between stiffness and tension, the mathematical model itself may be extended further to a problem of instability formation which is closer to hydraulic fracture. 81 Bibliography [1] Instability of Flows, volume 3 of Course of Theoretical Physics. Pergamon Press, Oxford; New York, 1989, c1977. [2] A. P. Bunger, E. Detournay, and D. I. Garagash. Toughness-dominated hydraulic fracture with leak-off. International Journal of Fracture. Reviews of Modern Physics, 134(2):175–190, 2005. [3] W. K. Chan and C. Yang. Surface-tension-driven liquid-liquid displacement in a capillary. Journal of Micromechanics and Microengineering, 15:1722–1728, 2005. [4] J. A. Diez, L. Kondic, and A. Bertozzi. Global models for moving contact lines. Physical Review E, 63, 200. [5] E. B. Dussan. On spreading of liquids on solid surfaces: Static and dynamic contact lines. Annual Review of Fluid Mechanics, 11:371–400, 1979. [6] R. Huang and Z. Suo. Wrinkilng of a compresssed elastic film on a viscous layer. Journal of Applied Physics, 91:1135–1142, 2002. [7] L. Kondic. Instability in Gravity Driven Flow of Thin Fluid Films. SIAM Review, 45(1):95–115, 2003. [8] L. Kondic and J. Diez. Pattern formation in the flow of thin films down an incline: Constant flux configuration. SIAM Review, 13(11):3168– 3185, 2001. [9] A. A. Lacey, J. R. Ockendon, and A. B. Tayler. “Waiting-time” solutions of a nonlinear diffusion equation. SIAM, 42(6), 1982. [10] T. G. Myers. Thin films with high surface tension. SIAM Review, 40(3):441–462, 1998. 82 [11] A. Oron, S. H. Davis, and S. G. Bankoff. Long-scale evolution of thin liquid films. Reviews of Modern Physics, 69:931980, 1997. [12] S. M. Roper and J. R. Lister. Buoyancy-driven crack propagation: the limit of large fracture toughness. Journal of Fluid Mechanics, 580:359– 380, 2007. [13] A. C. Slim, N. J. Balmforth, R. V. Craster, and J. C. Miller. Surface wrinkling of a channelized flow. Proceedings of Royal Society A, 465:123–142, 2009. [14] D. A. Spence, P.W. Sharp, and D.L. Turcotte. Buoyancy-driven crack propagation: a mechanism for magma migration. Journal of Fluid Mechanics, 174:135–153, 1987. [15] E. O. Tuck and L.W. Schwartz. A numerical and asymptotic study of some third-order ordinary differential equations relevant to draining and coating flows. SIAM Review, 32:453–469, 1990. [16] I. Veretennikov, A. Indeikina, and H.C. Chang. Front dynamics and fingering of a driven contact line. Journal of Fluid Mechanics, 373:81– 110, 1998. 83 Appendix A Numerical Details A.1 Horizontal Flow: Initial Condition The initial conditions of the fluid bump are modeled by the equation 3.8 ρg P h+ . (A.1) B B This is a linear equation and it can be solved analytically, giving h(x) as a function of exponentials and trigonometric functions. hxxxx = P h = A exp(βx) + B exp(−βx) + C cos(βx) + D sin(βx) − , ρg r β= 4 ρg B However, the next step requires imposing boundary conditions at a boundary that is unknown (contact point at θ is an unknown parameter in the formulation). This results in a set of equations which certainly cannot be solved analytically. Thus, either the original equation with the imposed boundary conditions or a system of seven nonlinear equations has to be solved numerically. Here, the preference was given to the former and the equation was solved using MATLAB built-in function bvp4c. The boundary conditions were chosen as described in Chapter 3. The system of differential equations was written as: y1 (x) y2 (x) y3 (x) y4 (x) = = = = h(x), h0 (x), h00 (x), h000 (x), Z r(θ) y5 (x) = h2 dx. (A.2) 0 84 A.2. Numerical Solution of the PDE Rewrite A.2 as a system of differential equations: y10 (x) = y2 (x), y20 (x) = y3 (x), y30 (x) = y4 (x), y40 (x) 4 (A.3) = [r(θ)] ρg P y1 (x) + B B , y50 (x) = y12 . Here, x has been re-scaled so the unknown parameter r(θ) is equal to unity when it appears in the boundary conditions. Next, MATLAB solves for h(x) on the interval (0, 1) and for the values of P and θ. The bvp4c solver also requires an initial condition (initial guess). The most effective way to select an initial condition is to solve the equation once with a fixed reasonable θ and use the obtained solution as the initial guess. The reason for this is that undefined boundary is a numerically difficult problem to deal with and if the initial guess is not good enough, the software may run into convergence issues. Pre-solving it once stabilizes the method. A.2 Numerical Solution of the PDE The equation to be solved is Equation 3.5, restated below ht = h3 hxxxxx + h3 hx x . (A.4) The Partial Differential Equation 3.5 is solved via Finite Differences. The accuracy of the method is determined by higher order derivatives since their numerical approximation requires the largest number of sample points. A.2.1 Time Discretization The time derivative in the model is first order and can be expressed as a first order derivative approximation. Since explicit schemes tend to be very sensitive to larger time steps, especially for such a high order problem, an implicit scheme is considered. Backwards Euler gives hk+1 − hk = A(hk+1 ), ∆t k = 1, ..., Nt − 1, (A.5) 85 A.2. Numerical Solution of the PDE where A(h) is the right hand side evaluated at h. Since A(h) is not a linear function, equation (A.5) cannot be re-arranged to solve for hk+1 in terms of hk and an iterative method has to be employed instead. Take the current approximation of hk+1 at the nth step to be hn so, hn+1 = hn + δhn . The goal is to solve for the correction term δhn by using Newton’s Approximation on A(h), hn + δhn − ∆t[A(hn ) + A0 (hn )δhn ] = hk . [I − ∆tA0 (hn )]δhn = rn where rn = hk − [hn − ∆tA(hn )]. (A.6) Using this method, each iteration solves for δhn and then hn +δhn is used as the guess in the next iteration. This is repeated until hn + δhn makes the residual term rn sufficiently small, so hk+1 = hn+1 . A.2.2 Space Discretization In order to apply the iteration method described, A(h) and A0 (h) have to be evaluated numerically at the mesh points. Consider the discretization of the domain 0 ≤ x ≤ L as xi = i∆x with i = 0...Nx and ∆x = L/Nx . Then the time derivative at each point is of the form dh(xi , t) = Ai , i = 1, ..., Nx − 1, (A.7) dt where Ai = A(h(xi , t)) is the discretization of the right hand side of equation (A.4). The left hand side has been described earlier, so this sections focuses on approximating the Ai term, A(h(xi , t)) = [h(xi , t)]3 [h(xi , t)]xxxxx + [h(xi , t)]3 [h(xi , t)]x x , or, to simplify the notation, 86 A.2. Numerical Solution of the PDE Ai = [hi ]3 [hi ]xxxxx + [hi ]3 [hi ]x x . The derivatives can be approximated via Finite Differences. The standard approaches are either backward difference and forward difference methods, both of order ∆x accurate, or the central differences, which is of order ∆x2 accurate. While one might be tempted to use the central differences because of higher accuracy, it turns out that it is not the most effective way to approximate a higher order derivative. Finite differences require a very large stencil (number of mesh points used in approximating the derivative) and for a sixth order derivative central differences would require six mesh points on each side, that is, a computational stencil of thirteen. The first order methods require only one extra point on either side (depending if it is forward or backward), so alternating between the two, the computational stencil size is only seven. Define the differences as hi − hi−1 hi+1 − hi , hi,x̄ = . ∆x ∆x Then the numerical discretization of Ai is hi,x = (A.8) Ai = [hi ]3 [hi ]x̄xx̄xx̄ + [hi ]3 [hi ]x̄ x . Notice that the h3i term is only differentiated once, so to maximize the symmetry of the problem we approximate it as the average of two consecutive h3 +h3 points. That is, h3i = i 2 i−1 , so that 3 hi + h3i−1 h3i + h3i−1 [hi ]x̄xx̄xx̄ + [hi ]x̄ . (A.9) Ai = 2 2 x The last task in the numerical approximation is to consider the derivative of Ai . The Jacobian matrix is then of the form, ∂A ∂A1 ∂A2 · · · ∂hN1x ∂h1 ∂h1 ∂A1 ∂A ∂A2 ∂h2 · · · ∂hN2x ∂h2 . .. .. .. . . . . . . ∂A1 ∂hNx ∂A2 ∂hNx ··· ∂ANx ∂hNx The derivatives can be computed from equation (A.9) with slight adjustments coming from the boundary conditions. 87 A.2. Numerical Solution of the PDE A.2.3 Numerical Boundary Conditions The original PDE boundary conditions given in equation (3.6) can be restated numerically, treating h(0, t) = h1 and h(L, t) = hNx h1 = δ, h2 − h0 = 0, 2∆x h2 − 2h1 + h0 = 0, 2∆x hNx = δ, hNx +1 − hNx −1 = 0, 2∆x hNx +1 − 2hNx + hNx −1 = 0. 2∆x (A.10) From equation (A.10) we can solve for the values of “ghost” points, h0 and hNx +1 , in terms of real mesh points and obtain a closed form expression for the neighbouring points, h2 and hNx +1 . These values are used to solve for Ai ’s close to boundaries. That is, everything before A4 and after ANx −2 . A concern one might have is the validity of using a pre-wetted thickness δ. The main concern comes from the fact that different δ’s would have a different effect on the flow and thus, resulting in different shape. However, this is not the case as Figure A.1 demonstrates. Varying δ by a factor of 10 still does not change the shape and the only difference between simulations is that at larger pre-wetted thicknesses more volume is added to the flow. Thus, as long as δ is chosen sufficiently small, the simulation will not have any dependence on that parameter. 88 A.2. Numerical Solution of the PDE 7 t0 = 0 6 5 Height t1 > t0 4 t2 > t1 t3 > t2 3 2 δ = 0.1 δ = 0.5 δ =0.01 1 0 0 5 10 15 20 Width 25 30 35 40 Figure A.1: The simulation of the horizontal flow with three δ values varying from 0.01 to 0.1. Each simulation run is graphed at four different times. A.2.4 Extending Numerics for Solving the PDE The differences between PDEs 3.5 and 4.9 is that the latter has an extra term present and a non-dimensional constant D. The structure of the model does not change, hence there are no changes to the time discretization A.6. Instead, a different A has to be used: h3i + h3i−1 3 3 [hi ]x̄ + [hi ] [hi ]x̄ . Ai = [hi ] [hi ]x̄xx̄xx̄ + D 2 x Also, A0 (hn ) and the boundary conditions are adjusted accordingly. The non-dimensional constant D adds an additional degree of freedom. During the numerical analysis this constant is set to a reasonable quantity with an order of magnitude of one. In Chapter 5 the simulations are repeated with the appropriate value, as determined by the physical quantities. 89 A.3. Inclined Plane: Solution of the ODE 4.11 A.3 Inclined Plane: Solution of the ODE 4.11 This non-linear ODE cannot be solved analytically. However, there is enough information known about the behaviour of the solutions that this boundary value problem can be solved numerically via the MATLAB built-in solver bvp4c, provided the ODE is kept in one dimension. The governing equation is U f + f 3 fξξξξξ + Df 3 fξ − f 3 = c, (A.11) with the boundary conditions defining the constants U = 1 + δ + δ2, c = δ + δ2. The boundary conditions also require that the left hand side of the solution remains at a constant level of unity as the fluid flows to the right and drops off to the pre-wetted thickness δ as the solution reaches the right boundary. Observe that the constant solutions f = δ and f = 1 both satisfy A.11, but do not satisfy the boundary conditions. Thus, we need to find a solution which travels from f = 1 at the left boundary to f = δ at the right boundary. It is also expected that the derivative of f should vanish as it approaches the required constant level. A natural way to model this is to consider an exponential perturbation, which would grow when f is close to one (positive eigenvalues) and diminish when f is close to δ (negative eigenvalues). In order to find these eigenvalues, consider a homogeneous linearization of equation A.11 about a constant, µ, with a small perturbation ∆ = eλξ : 3(U − 1) − 2U µ = 0. ∆ξξξξξ + D∆ξ + ∆ µ4 Then at the left boundary, when f = 1, we must take the positive roots of the characteristic equation ∆ξξξξξ + D∆ξ + (U − 3)∆ = 0. Similarly, the right boundary conditions need the negative roots of the characteristic equation 90 A.3. Inclined Plane: Solution of the ODE 4.11 U − 3δ 2 = 0. δ3 Then the solution at the left boundary has the form ∆ξξξξξ + D∆ξ + ∆ f =1+ n X aj e λ j ξ . j=1 This allows us to solve for the unknown constants aj in terms of f and its derivatives and consequently to re-express the higher order derivatives in terms of the lower order derivatives. Enforcing the higher order derivatives to vanish stabilizes the problem while still imposing the required behaviour on the solution. A.3.1 Growth Rate Eigenvalue Solution The solution of the eigenvalue problem 4.25 uses the same methodology as the solution of the ODE. The growth and decay conditions are enforced at the boundaries. The eigenfunction φ behaves like a derivative of h, so we now linearize about 0 at both boundaries, giving the new characteristic equation 0= h30 λ6 − 3q 2 h30 λ4 4 + (3q + D)h30 λ2 + 3c − 2U h0 λ + q 2 h30 (q 4 − D) − σ. The rest of the procedure is almost identical that already described. 91 Appendix B Experimental Set Up Materials: Experiments are done at SATP. The experimental fluid is corn syrup (viscosity ∼ 8.5 ± 0.5 Pa · s at room temperature). The elastic latex sheet (thicknesses 0.15 ± 0.02 mm and 0.20 ± 0.02 mm, Young’s modulus 2 ± 0.4 GPa) is purchased from McMaster Carr. The results presented here are from experiments with latex of thickness 0.20 mm. The bending stiffness B = 0.00533N · m is determined by integrating over the thickness of the the film (in SI units). Z 0.0020 2 · 10−9 z 2 dz B= 0 Corn syrup is dyed with black ink at 1 mL of ink per 25 mL of fluid. Procedure: The fluid (60-100 mL) is placed on a clear plexiglas (39 by 78 cm) baseboard and the latex sheet is taped on top. The exact volume is determined by how large of a latex sheet is used. To ensure no air bubbles lie between the latex and the fluid, the liquid must completely cover the surface area under the sheet. Thus, thin pre-wetting is present in the experiment as well as the numerical simulations. A light source is placed below the plexiglas base and is diffused with a white screen (see Figure B.1). The images of the flow are taken from the top and the thickness of the fluid at a point is measured from the light intensity transmitted through the apparatus at that point. In the case of flow on a perfectly horizontal baseboard, the fluid is collected in the centre and then allowed to spread out freely. 92 Appendix B. Experimental Set Up Camera 78 cm Fluid 39 cm 24 cm White Screen to Diffuse Light 46cm Light Source Figure B.1: Experimental set up Calibration: In order to measure the height of the wave front, the thickness has to be calibrated. This is done by using a wedge of known height (from 0 mm to 8.53 mm) filled with corn syrup and then calculating the calibration curve from the measurements along this wedge. When the direction of the flow is parallel to the light source, the wedge can be left on the side of the baseboard so that the thickness of the spreading fluid can be determined at any point in time. 93 Appendix C MATLAB Code Here we present the core code that was used to solve the partial differential equations 3.5 and 4.10 via Newton’s iteration. The solution makes calls to functions Ahn6 and Ahnpa6 which compute the righthand side of the equations and its derivative at each time interval. These computations differ slightly depending on the inclination. The call made to movie movieContour, which retrieves the first frame of the experimental video recording. %SolutionPDE %Retrieving leftP = 54; [height xh] xc = 24.96; experimental initial conditions rightP = 75; %boundaries of the fluid position = movieContour; hc = 5.9; %characteristic parameters %program relevant parameters N = 200; c=(xh(rightP)-xh(leftP))*4; a=0; b = 0.01; b1 = b; d = 1; x = (c-a)/N; X = a:x:c; n = length(X); maxit = 100; tol = 10^(-10); ep =sqrt(eps); I = eye(n); height = height/hc; %Polynomial fit of the initial condition and imposing symmetry hic = spline(xh(leftP:rightP)-xh(leftP), ... height(leftP:rightP), X(1:(end+3)/4)); b = hic(end); b1 = b; h = b*ones(1,length(X)); 94 Appendix C. MATLAB Code h((end-1)/4:3*(end+3)/4-1) = [hic(end:-1:1) hic(1) hic]; h = h’; X = X/xc; % % Numerical boundary conditions for incline plane % h = 1/2*(1+b)-1/2*(1-b)*tanh(2X - 10)-(1-b)*exp(-(X.^2)/10); % h = h’; %time related parameters t = 0; dt = 0.1; tmax = 2000; tplot = 5; nplots = round(tmax/tplot); plotgap = round(tplot/dt); dt = tplot/plotgap; conv = 0; %counter for non-convergent loops total = 0; %counter for average # of iterations %code if producing a video of the simulation % plot(X,h); % axis([0 20 0 2]); % Hmovie(1)= getframe; plot(X,h); hold all; drawnow; %Solving the PDE for i = 1:nplots for j = 1:plotgap %two loops are used so the solution can %be plotted only at selected intervals t = t+dt; z = (i-1)*plotgap + j; hk = h; A = Ahn6(h,x,n,d); R = dt*(A); %newton’s method to solve for h for each timestep for it = 1:maxit 95 Appendix C. MATLAB Code J = (I - dt*Ahnpa6(h,x,n,d)); dH = J\R; h = h + dH; A = Ahn6(h,x,n,d); R = hk - h + dt*A; RS(it) = norm(R,2); if RS(it) < tol*norm(h,2); break; end end if it >= maxit, msg=[’WARNING: timestep ’, num2str(z), ... ’ did not converge’] conv = conv +1; end total = total +it; % Determining parameters: front position and maximum height Hbump(z) = max(h); %interpolate the peak position, so it is not rounded %to nearest mesh point sr = find(hr <= b,1); pp = polyfit(X(sr-2:sr),[hr(sr-2) hr(sr-1) hr(sr)],2); Xr(z) = -pp(2)/2/pp(1) - sqrt(pp(2)^2 - ... 4*pp(1)*(pp(3)-b))/2/pp(1); % impose boundary conditions h(end) = b; h(1) = b1; h(end-1) = b; h(2) = b1; end Hploted(:,i) = h; 96 Appendix C. MATLAB Code plot(X,h); hold all; drawnow; %code if producing a video of the simulation %plot(X,h); %axis([0 20 0 2]); %Hmovie(i+1)= getframe; end %movie(F); % displays information if all the iterations converged % and how quickly msg = [num2str(conv),’ out of ’, num2str(i*j),... ’ loops did not converge’] msg = [’Average # of iterations ’, num2str(total/(i*j))] 97 Appendix C. MATLAB Code The following Ahn6 and Ahnpa6 functions were used when solving equations 4.10. The changes required to modify this code for 3.5 are indicated as comments in the code. % % % % Accepts nx1 vector h, a real x (equidistant mesh size) and integer n (# of mesh points) non-dimensional constant d Returns nx1 vector corresponding to A(h) function A = Ahn6(h,x,n,d) %matrices approximating the derivatives s0 = ones(n,1); s1 = ones(n-1,1); s2 = ones(n-2,1); s3 = ones(n-3,1); D41 = sparse(diag(s3,3) + diag(-5*s2,2) + diag(10*s1,1)... + diag(-10*s0) + diag(5*s1,-1) + diag(-s2,-2)); D42 = sparse(diag(s3,3) + diag(-6*s2,2) + diag(15*s1,1)... + diag(-20*s0) + diag(15*s1,-1) + diag(-6*s2,-2)... + diag(s3,-3)); D43 = sparse(diag(s2,2) + diag(-5*s1,1) + diag(10*s0)... + diag(-10*s1,-1) + diag(5*s2,-2) + diag(-s3,-3)); D21 = sparse(diag(s1,1) + diag(-s0)); D22 = sparse(diag(s1,1) + diag(-2*s0) + diag(s1,-1) ); D23 = sparse(diag(s0) + diag(-s1,-1)); D1 = sparse(diag(s1,1) + diag(-s1,-1)); Du = sparse(diag(s1,1)); Dd = sparse(diag(s1,-1)); %computation of A with above matrices h3 = h.^3; A4 = ((Du*h3).*(D41*h) + (h3).*(D42*h)... - (Dd*h3).*(D43*h))/(2*x^6); A2 = d*((Du*h3).*(D21*h) + (h3).*(D22*h)... - (Dd*h3).*(D23*h))/(2*x^2); A1 = -(D1*h3)/(2*x); %A1 = 0 %for solving the PDE with 0 inclination 98 Appendix C. MATLAB Code A = A4 + A2 + A1; A(1) = 0; A(2) = 0; A(end) = 0; A(end-1) = 0; A(3) = A(3) + ((h3(3)+h3(2) )*h(2))/(2*x^6); A(end-2) = A(end-2) + ((h3(end-1)+h3(end-2))*h(end-1))/(2*x^6); 99 Appendix C. MATLAB Code % % % % Accepts nx1 vector h, a real x (equidistant mesh size) and integer n (# of mesh points) non-dimensional constant d Returns nx1 vector corresponding to A’(h) function A = Ahnpa6(h,x,n,d) %matrices approximating the derivatives s0 = ones(n,1); s1 = ones(n-1,1); s2 = ones(n-2,1); s3 = ones(n-3,1); h2 = h.^2; h3 = h.^3; D41 = sparse(diag(s3,3) + diag(-5*s2,2) + diag(10*s1,1)... + diag(-10*s0) + diag(5*s1,-1) + diag(-s2,-2)); D42 = sparse(diag(s3,3) + diag(-6*s2,2) + diag(15*s1,1)... + diag(-20*s0) + diag(15*s1,-1) + diag(-6*s2,-2)... + diag(s3,-3)); D43 = sparse(diag(s2,2) + diag(-5*s1,1) + diag(10*s0)... + diag(-10*s1,-1) + diag(5*s2,-2) + diag(-s3,-3)); D2 = sparse(diag(s1,1) + diag(-2*s0) + diag(s1,-1)); D2u = sparse(diag(s1,1) + diag(-s0)); D2d = sparse(diag(s0) + diag(-s1,-1)); D1 = sparse(diag(s1,1) + diag(-s1,-1)); Du = sparse(diag(s1,1)); Dd = sparse(diag(s1,-1)); h3u = Du*h3; h2u = Du*h2; h3d = Dd*h3; h2d = Dd*h2; A0 = ( -10*h3u - 20*h3 - 10*h3d + 3*h2.*(D42*h) )/(2*x^6)... + d*(3*h2.*(D2*h)-h3d-h3-h3u)/(2*x^2); Ap1 = ( 10*h3u + 15*h3 + 5*h3d + 3*h2u.*(D41*h) )/(2*x^6)... + d*(h3u+h3+3*h2u.*(D2u*h))/(2*x^2) ... - 3*h2u/(2*x); %this last term is removed when %solving the PDE with 0 inclination Am1 = ( 5*h3u + 15*h3 + 10*h3d - 3*h2d.*(D43*h) )/(2*x^6)... + d*(h3d+h3-3*h2d.*(D2d*h))/(2*x^2) ... 100 Appendix C. MATLAB Code + 3*h2d/(2*x); %this last term is removed when %solving the PDE with 0 inclination Ap2 = (-5*h3u - 6*h3 - h3d)/(2*x^6); Am2 = (-h3u - 6*h3 - 5*h3d)/(2*x^6); Ap3 = ( h3u + h3 )/(2*x^6); Am3 = ( h3 + h3d )/(2*x^6); A = sparse(diag(Ap3(1:end-3),3) + diag(Ap2(1:end-2),2) ... + diag(Ap1(1:end-1),1) + diag(A0) + diag(Am1(2:end),-1) ... + diag(Am2(3:end),-2) + diag(Am3(4:end),-3)); %Boundary Conditions A(1,:) = 0; A(2,:) = 0; A(end,:) = 0; A(end-1,:) = 0; A(3,2) = A(3,2) + (4*h3(2)+h3(3) )/(2*x^6); A(3,3) = A(3,3) + 3*h2(3)*h(2)/(2*x^6); A(end-2,end-1) = A(end-2,end-1)+(h3(end-2)+4*h3(end-1))/(2*x^6); A(end-2,end-2) = A(end-2,end-2)+3*h2(end-2)*h(end-1)/(2*x^6); % % Boundary Conditions for solving the PDE with 0 inclination % A(1,1) = A(1,1) + (- 10*h3(2) + 3*h3(1)*(15*h(2) - 6*h(3) ... + h(4)) )/(2*x^6) + (3*(h(1))^2*h(2)- h3(2) )/(2*x^2); % A(1,2) = A(1,2) + (80*h3(2) - 7*h(3) -10*h(1) + h(4) ... +15*h3(1))/(2*x^6) + (h3(1) - 3*(h(2))^2*h(1) + 4*h3(2))/(2*x^2); % A(1,3) = A(1,3) + ( -7*h3(2)-6*h3(1) )/(2*x^6); % A(1,4) = A(1,4) + ( h3(2) + h3(1) )/(2*x^6); % % A(2,1) = A(2,1) + (-15*(h(1))^2*h(2)+3*(h(1))^2*h(3))/(2*x^6); % A(2,2) = A(2,2) + (-h3(3)-24*h3(2)+3*(h(2))^2*h(3) ... % - 5*h3(1))/(2*x^6); % A(2,1) = A(2,1) + ( h3(2)+h3(1)-3*(h(3))^2*h(2))/(2*x^6); % % A(end,:) = 0; A(end-1,:) = 0; A(end-2,:) = 0; % % A(3,2) = A(3,2) + (4*h3(2)+h3(3) )/(2*x^6); % A(3,3) = A(3,3) + 3*h2(3)*h(2)/(2*x^6); 101
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Viscous fluid instabilities under an elastic sheet
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Viscous fluid instabilities under an elastic sheet Khomenko, Maria 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Viscous fluid instabilities under an elastic sheet |
Creator |
Khomenko, Maria |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | This thesis considers the flow of thin fluid film between an elastic sheet and a rigid plane. We derive a mathematical model for the flow from the Navier-Stokes equations using the lubrication approximation and develop numerical and similarity solutions to this problem. An experimental apparatus was developed to investigate this phenomenon, and the results of the mathematical model were compared with experimental data. Chapter 3 examines the evolution of a fixed fluid volume under gravitational forces on a horizontal plane. The evolution of the fluid mass profile and the progression of the fluid front are determined from the numerical solutions, as well as experimentally. The favourable comparison between the numerical solutions and the experimental results establishes the validity of the model. Chapters 4-5 considers the evolution of a thin fluid flow under an elastic on an inclined plane. We establish a traveling wave solution for this flow. A linear stability analysis yields the criterion for the existence of unstable modes and establishes the growth rate and wavelength of the most unstable mode. Instability is promoted by increasing the inclination of the plane. For low angles, the numerical and experimental growth rates were in good agreement, while the wavelengths were experimentally of the same order and numerically computed wavelengths had little variation. The long term behaviour of the fluid front is studied analytically via a similarity solution in Chapter 6. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0069934 |
URI | http://hdl.handle.net/2429/23813 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_khomenko_maria.pdf [ 1.7MB ]
- Metadata
- JSON: 24-1.0069934.json
- JSON-LD: 24-1.0069934-ld.json
- RDF/XML (Pretty): 24-1.0069934-rdf.xml
- RDF/JSON: 24-1.0069934-rdf.json
- Turtle: 24-1.0069934-turtle.txt
- N-Triples: 24-1.0069934-rdf-ntriples.txt
- Original Record: 24-1.0069934-source.json
- Full Text
- 24-1.0069934-fulltext.txt
- Citation
- 24-1.0069934.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0069934/manifest