Simulations of Falling, Flexible Paper

Presented at the SC99 Conference 11/13-19/1999

Comer Duncan and Guangnian Zhai
Department of Physics and Astronomy
Bowling Green State University
Bowling Green, OH 43403


The computations reported here were performed on the OSC Cray T-94 and we thank the Center for its computational support. The visualizations were performed using AVS Express, available at OSC. We especially thank Ken Flurchick for sharing his expertise on AVS and visualization in general as well as making the mpeg movies of our simulations.


This report is work in progress on our simulations of the dynamics of falling, flexible paper. The results discussed herein are thus preliminary and represent some of our results to date. Further modifications will be made to the methods and modes of analysis. The final results are to be reported in a forthcoming paper.

I. Description of Project Problem and Previous Approaches

The understanding of the dynamics of falling objects is a long standing area of interest. When the falling object is not rigid, such as in the case of a falling piece of paper, the issues of its dynamics and the relation of the dynamics to that of the fluid in which it falls has a long and interesting history going back to works of Maxwell [1], Helmholtz, Kelvin, and others [2,3]. Due to the fundamental coupling between the immersed object and the fluid in which it falls, understanding dominant aspects of the motion of the falling paper has proven to be difficult. In the last two years, there have been new experiments performed to explore and classify the types of motion of falling strips of various materials. The experimental work of Belmonte et al [4] and that of Mahadaven et al [5] have resulted in advances in understanding the dependence of the types of paper motion on the initial conditions and on the Froude number. It has been found that the motion exhibits a transition from a fluttering to a tumbling behavior at a critical value of the Froude number. Nevertheless, many gaps remain in fundamentally understanding how the interplay of the motion of the flexible paper and the incompressible fluid in which it falls determine the complete motion.

Theoretically understanding the motion of the paper has been the subject of a number of works over the past few years. The work of Tanabe and Kaneko [6] proposed a phenomenological model in which the paper was subject to lift, friction, and gravity. They treated the paper as a rigid object and developed a dynamical system including the forces of lift and a phenomenological model for drag. The resulting system of ordinary differential equations was numerically integrated in various limits. The system was found to exhibit several different classes of behaviors including periodic rotation, chaotic fluttering, periodic fluttering, and a simple perpendicular fall. The theoretical models put forward by Belmonte et al [4] in conjunction with their experiments were also built upon a phenomenological model for the drag on the paper due to the fluid. Both of these works did not attempt to model the paper and the fluid in which it is immersed in a self-consistent manner, electing to introduce phenomenological elements presumably characterizing the essential aspects of the fluid dynamics. In addition, the experiments of Belmonte et al also visualized the motion of the fluid in which the strips fall, giving ample evidence for quite complex fluid dynamics at and around the falling strips. It is notable that both experiments utilized rather stiff strips of various materials and that the fluid in which the strips fell were liquids. By restricting to stiff materials, the extra complexities which would be present due to flexibility were minimized, perhaps allowing a more controlled examination of the dependence of the motion on Froude number. However, it is clear that 'real' paper is flexible and that the motion of flexible materials will possess extra complexity. We feel that the inclusion of flexibility is an important step toward the goal of a completely characterizing the dynamics of both the paper and the fluid in which it falls.

However, if the paper is allowed to be flexible, the numerical simulation's complexity of the paper's dynamics will increase. This is due to several factors. First, when the paper is flexible the paper itself will exert forces back on the fluid which are perhaps different depending on where along the paper one examines the force. The forces at the ends of the paper are different than at interior points. Second, the effect of flexible paper on the surrounding fluid is then potentially quite complex, with the likely induced vorticity of the fluid and its back reaction on the paper. As the paper flexes and falls, the vortices produced will be shed and the back reaction on the flexible paper will influence its motion. This makes the use of phenomenological models of drag and lift of doubtful utility. In fact the approach we propose below would not include any assumed phenomenological forces; those effects would emerge from a proper description.

With these considerations in mind, we feel that what is needed is a totally self-consistent model in which both the flexible paper and the fluid in which it is immersed are treated on an equal dynamical footing. Such a need presents us with a challenge to treat the paper dynamics and the fluid dynamics without any ad hoc assumptions about lift or drag.

The present preliminary research report, conducted using OSC facilities, gives an account of our results to date of such a coupled treatment of paper and fluid dynamics, which can capture the relevant types of motion using as close to a first principles treatment as possible. In addition, we report preliminary visualization results of the induced motion of fluid. We can track the development and subsequent evolution of vortices and their shedding for several different systems. A long range goal is to assess the drag for this strongly coupled system.

With the successful completion of the project, we expect to have the first computational treatment of the falling, flexible paper problem in which we can determine the effect of the fluid's dynamics on the motions of falling paper. After a hundred years, we are now in a position to investigate the complete dynamics of this interesting and complex system and make contact with recent experiments. It is hoped that such work will contribute some additional elements to fundamentally understand the system and perhaps suggest further experiments to probe the interesting dynamics of this strongly coupled paper-fluid system.

II. Computational Methods

The approach we will use to perform these simulations stems from work done by Charles Peskin and David McQueen [7-10] at Courant Institute of Mathematical Sciences in their development of a computational model of the heart. Peskin [7] developed a method called the Immersed Boundary Method which treats the heart as an immersed boundary in the fluid and allows a close coupling between the fluid and the flexible boundary that the heart constitutes. A set of particles model the fibers which make up the heart. The fibers interact with the blood and exert forces on the blood. The blood exerts forces on the flexible immersed boundary fibers inducing the walls to bend. The net system is quite complex. Over the years, Peskin and McQueen have made major strides in such modeling so that by now the Immersed Boundary Method is well developed and has been shown to be robust, providing a sophisticated three dimensional computational model which compares quite well to the properties of real mammalian hearts.

It is clear that this method can be applied to a variety of systems in which the immersed boundary is part and parcel of the appropriate description. The falling paper problem is clearly within the domain encompassed by the Immersed Boundary Method and thus we have adopted it as the primary computational tool in our approach.

The present computations initially have used our own two dimensional immersed boundary method code. We have written the code using methods described in publications of Peskin and McQueen [7-10]. The code includes a Navier-Stokes solver and assumes periodic boundary conditions along each rectangular boundary. The current code uses Chorin's method [11] and has been used by other workers [12]. We have performed simulations using this code on workstations at Bowling Green State University and on the OSC Cray T-94.

The fluid's dynamical influence on the paper is a major aspect of this work. As such, the flow's characterization in the system is of utmost importance and the output will focus on these runs. We have done some preliminary work on visualizing the motion of the particles constituting the paper. We use AVS to perform the visualizations. Since the current project is in two dimensions, AVS is a quite useful tool to add ways to view the evolving flow as we consider different effective Froude numbers and study the relations between the apparent dependence of the complexity of the paper's and of the fluid's motions.

The coupled system of particles and fluid may be briefly characterized as follows. For a more complete description, one should consult several references of work by Peskin and McQueen, e.g. [10]. The fluid is incompressible. The fluid equations are thus taken to be the Navier-Stokes equations for the velocity field with the constraint that the velocity field is divergence-free. The Navier-Stokes equations for the velocity field is supplemented with a forcing term on its right hand side. This is the force on a given fluid cell due to the immersed boundary points. The force is mostly zero except where the particles are adjacent to the given cell of interest. The force is modeled by an approximate delta function so that the force on a given cell is effectively due to only those particles within a restricted distance from the given cell.

The particles which exist in the fluid in general do not reside at fluid cell positions. They are strung throughout the fluid. The dynamics of the particles are obtained by enforcing the no-slip boundary condition on the fluid velocity at the locations of the particles. The particle positions are Lagrangian points which move with the particles and the fluid velocity are interpolated to the particle positions. Note that in the Immersed Boundary Method the fluid is given an Eulerian description while the particles are given a Lagrangian description. The method handles the need to go back and forth between these modes.

In the form in which our code is written, the Chorin solver is first order accurate. The particle positions are advanced in time using the appropriately interpolated fluid velocity to the particle positions and a forward Euler method is used to perform the particle position updates.

III. Immersed Boundary Algorithm

The algorithm used in our present two-dimensional code may be described in the following steps for a single cycle through the code:

  1. compute the force on the immersed boundary points due to interparticle interactions
    force on kth particle fk = function(Xk,Xk+1, Xk-1) due to stretching and bending
  2. spread the boundary force to the fluid lattice points
    The force on the i,j lattice point
    F(i,j) = Sum over k=1,Nparticles (fk dirac-delta (x(i,j) - Xk(i,j))
  3. solve the Navier-Stokes equations with Force(i,j) as a source term. This gives new values of velocity components ux(i,j) and uy(i,j)
  4. interpolate the resulting velocity on the fluid lattice points to the immersed boundary points
    Uk = Sum over all i,j (u(i,j) dirac-delta (x(i,j) - Xk(i,j))
  5. update the boundary points to the new positions
    Xknew = Xkold + dt * Uk

Thus each time step there are two activities involving the delta function and one call to the Navier-Stokes solver. At the end of the process, the particle positions are updated.

IV. Computational Parameter Choices

Since the present code to simulate the falling paper is two dimensional we are actually simulating a falling stick which can be flexible. The experiments conducted have been constrained so that the motion is substantially two dimensional, so we feel that this first approach using the two dimensional code is appropriate.

The paper is approximated by a set of particles. In the simulations we have settled on 256 particles. These particles interact with each other with a spring-like force which tends to keep the particles at most a distance h/2 apart, where h has been set to the grid spacing of the computational mesh associated with the fluid. This stretching force is characterized by a Hooke's law form with a spring constant associated with the stretching parameter S_s. Secondly, the paper's flexibility is modeled using a bending force which tends to keep the paper straight. The strength of the total force component is the bending parameter Sb. Finally, the force of gravity acts on the particles [but not the fluid]. These forces are computed as a function of the instantaneous positions of the particles. The end particles have forces calculated which take into account that the end particle has a neighbor on one but not both sides. The net result is a force between the particles which models the material. By choosing Ss and Sb in relation to the force of gravity, one can obtain a variety of effective material properties. We have explored part of this parameter space with our workstation simulations and have determined the ranges of Ss and Sb which give a stiff or flexible two-dimensional string of particles.

The fluid which interacts with the paper is modeled using Chorin's Navier-Stokes solver and is currently restricted to square domains. An important computational aspect of Chorin's method for the simulation of the fluid centers on using the Fast Fourier Transform. For the tests we report below, we have settled on 256 x 256 for the incompressible fluid simulation. We have performed test runs with the following sets of parameters.

Number of fluid cells: 256 x 256
Number of particles for flexible paper: 256
Size of Fluid Domain: 6 cm x 6 cm
Length of paper: 1 cm
Total Mass of paper: 0.7 g
Time step [not adequately small] dt : 1.0 e-04 - 1.0 e-06 sec
Stretching Coefficient: 1.0 e 08-1.0 e 10
Bending Coefficient: 1.0 e 08 - 1.0 e 10
Fluid Density: 1.0 g/cm3
Fluid Viscosity: 0.01 - 0.0015 g/cm sec
Number of time steps of duration dt: 3.5 e 04 - 3.2 e 06
For the simulations we report here, the T-94 times varied between 4 and 20 cpu hours.

V. Some Preliminary Results

We present some representative results for two computational cases:

  1. Falling, Stiff Paper at Initially Horizontal Orientation

    The initial state of the paper is at an angle of zero degrees relative to the Horizontal. There are 256 particles constituting the paper with the separation between particles much smaller than the maximum of h/2. The purpose of this run was to explore the dynamics of the paper when the stretching and bending spring constants were chosen quite large [both 1.0 e 10]. The force of gravity acts in the -y direction. The motion of the fluid and paper were computed using our immersed boundary code and we show, in Figures 1a and 1b, the AVS visualization of both the particles and the fluid. Fig. 1a shows the velocity vectors and illustrates vortex shedding and the influence of the fluid velocity field on the paper. We note the beginning of the production of vortices and their shedding, along with the coupled interaction of the rather stiff paper and the fluid. We plan to use such simulations to explore the drag and lift force dependence on the velocity field.

    Animations of data slice plus vector arrows and contour lines of the magnitude of the velocity

    flat2_line.jpg Starting at angle=0o
    Frame 2
    24.16 Kb
    flat2.jpg Starting at angle=0o
    Frame 2
    21.03 Kb
    flat40_line.jpg Starting at angle=0o
    Frame 40
    44.23 Kb
    flat40.jpg Starting at angle=0o
    Frame 40
    51.79 Kb
    flat80_line.jpg Starting at angle=0o
    Frame 80
    54.45 Kb
    flat80.jpg Starting at angle=0o
    Frame 80
    74.65 Kb

  2. Falling, Flexible Paper at an Initial Angle of 60 degrees

    In this simulation, we begin to explore the types of motion that paper can possess when we allow the paper to be flexible. This is shown in Figures 2a and 2b. The present simulations were done with a viscosity which is quite large, so the damping of the motion is an important factor in the long term motion. Here, too, the dynamics of the fluid velocity field is visualized along with the paper.

    Animations of data slice plus vector arrows and contour lines of the magnitude of the velocity

    angle2_line.jpg Starting at angle=60o
    Frame 2
    23.78 Kb
    angle2.jpg Starting at angle=60o
    Frame 2
    20.41 Kb
    angle40_line.jpg Starting at angle=60o
    Frame 40
    44.02 Kb
    angle40.jpg Starting at angle=60o
    Frame 40
    51.93 Kb
    anfle80_line.jpg Starting at angle=60o
    Frame 80
    59.63 Kb
    angle80.jpg Starting at angle=60o
    Frame 80
    69.45 Kb

  3. Falling, Flexible Paper at an Initial Angle of 60 degrees

    Further simulations in the flexible paper range exhibit quite complex behavior. Here we begin to explore the rich diversity of the dynamics of flexible paper and fluid motions when the paper is moderately flexible and compare with the stiff paper simulations. Animations 3a and 3b shows the dynamics of the velocity vector field for one such flexible paper case. The vortex shedding and backreaction of the velocity field on the falling paper is evidently less than in the stiff paper case. However, the concomitant increase in complexity of the motion of both fluid and paper necessitates higher resolution and a Navier-Stokes solver with as small numerical viscosity as feasible. These simulations are thus exploratory only and will be part of an extensive study of the many facets of the flexible paper dynamics likely to be exhibited.

    Animations 3a and 3b

  1. J. C. Maxwell, in The Scientific Papers of James Clerk Maxwell (Dover, New York,1890), pp. 115-118.
  2. M. J. Lighthill, Annual Reviews of Fluid Mechanics 1, 413 (1969).
  3. H. J. Lugt, Annual Reviews of Fluid Mechanics 15, 123 (1983).
  4. A. Belmonte, H. Eisenberg, and E. Moses, Physical Review Letters 81, 345 (1998).
  5. L. Mahadevan, W. S. Ryu, and A. D. T. Samuel, Physics of Fluids 11, 1 (1999).
  6. Y. Tanabe and K. Kaneko, Physical Review Letters 73, 1372 (1994).
  7. C. S. Peskin, Journal of Computational Physics 10, 252 (1972).
  8. C. S. Peskin, Journal of Computational Physics 25, 220 (1977).
  9. C. S. Peskin and D. M. McQueen, Journal of Computational Physics 81, 372 (1989).
  10. C. S. Peskin and D. M. McQueen, Symposia of the Society for Experimental Biology, 49, 265 (1995).
  11. A. J. Chorin, Mathematics of Computation 22, 745 (1968).
  12. L. J. Fauci, Journal of Computational Physics 86, 294 (1990).
  13. G. Zhai and C. Duncan, "Simulations of Falling Flexible Paper Using the Immersed Boundary Method", in preparation (1999).