Research

Flexible fibers in Stokes flow

Publications (most recent first):

A.-K. Tornberg and B. Engquist.
High order difference methods for wave propagation in discontinuous media.
Submitted.

A.-K. Tornberg, B. Engquist, B. Gustafsson and P. Wahlund.
A new type of boundary treatment for wave propagation.
Accepted for publication in BIT.

A.-K. Tornberg and K. Gustavsson.
A numerical method for simulations of rigid fiber suspensions.
Accepted for publication in J. Comput. Phys . Available online from the publisher.

A.-K. Tornberg.
Regularization techniques for singular source terms in differential equations.
In Laptev, A. (ed.), European Congress of Mathematics (ECM), Stockholm, Sweden, June 27-July 2, 2004. Zurich: European Mathematical Society (2005). (ISBN 3-03719-009-4, pages 477-499).

B. Engquist, A.-K. Tornberg and R. Tsai.
Discretization of Dirac Delta Functions in Level Set Methods.
J. Comput. Phys. 207:28-51, 2005.

A.-K. Tornberg.
Numerical Simulations of the Dynamics of Fiber Suspensions.
In Multiscale methods in Science and Engineering (Lecture Notes in Computational Science and Engineering) , editors B. Engquist et al., Springer-Verlag, 2005 (ISBN 3540253351).

M.J. Shelley and A.-K. Tornberg.
Microstructural Dynamics in Complex Fluids.
In the chapter on Mathematical Methods (editor M. Bazant) in Handbook of Materials Modeling .
Kluwer Academic Publishers, 2005 (ISBN 1402032870).

A.-K. Tornberg and B. Engquist.
Numerical approximations of singular source terms in differential equations
J. Comput. Phys. 200:462-488, 2004.

A.-K. Tornberg and M.J, Shelley.
Simulating the Dynamics and Interactions of Flexible Fibers in Stokes Flows
J. Comput. Phys. 196:8-40, 2004.

A.-K. Tornberg and B. Engquist.
Regularization Techniques for Numerical Approximation of PDEs with Singularities.
Journal of Scientific Computing . 19:527-552, 2003.

A.-K. Tornberg and B. Engquist.
The Segment Projection Method for Interface Tracking.
Commun. Pur. Appl. Math. , 56:47-79, 2003.

B. Engquist, O. Runborg and A.-K. Tornberg.
The Segment Projection Method for Geometrical Optics.
J. Comput. Phys. 178:373:390, 2002.

A.-K. Tornberg.
Multi-Dimensional Quadrature of Singular and Discontinuous Functions.
BIT . 42:644-669, 2002.

A.-K. Tornberg and B. Engquist.
Interface Tracking in Two-Phase Flows.
Multifield Problems, State of the Art , Springer Verlag, Berlin, 58-66, 2000.

A.-K. Tornberg.
Interface Tracking Methods with Application to Multiphase Flows. (This file is huge).
PhD Thesis , NADA, KTH, Stockholm, Sweden, 2000. ISBN 91-7170-558-9, TRITA-NA 0010.

A.-K. Tornberg and B. Engquist.
Finite Element Based Methods for Interfacial Flow Simulations.
Finite Element Based Methods for Interfacial Flow Simulations. Proceedings of ENUMATH99,
World Scientific Publ. 251-258, 2000.

A.-K. Tornberg and B. Engquist.
A Finite Element Based Level Set Method for Multiphase Flow Applications.
Computing and Visualization in Science, 3:93-101, 2000.

A.-K. Tornberg.
A Finite Element Based Level Set Method for Multiphase Flow Simulations.
Licentiate's Thesis , NADA, KTH, Stockholm, Sweden, 1998. ISBN 91-7170-317-9, TRITA-NA-9817.

A.-K. Tornberg, R.~W. Metcalfe, R.~Scott, and B.~Bagheri.
A Fluid Particle Simulation Method
In Computational Science for the 21st century (M-O. Bristeau et al., eds.). John Wiley & Sons, New York, pp 312-321 (1997).

A.-K. Tornberg, R.~W. Metcalfe, R.~Scott, and B.~Bagheri.
A Front-Tracking Method for Simulating Fluid Particle Motion using High-Order Finite Element Methods.
1997 ASME Fluids Engineering Division Summer Meeting FEDSM97-3493.

Interface tracking methods

My research is in the area of numerical methods for partial differential equations, and my main interest is moving boundary problems. Interfaces or internal boundaries are present in many different applications, such as high frequency wave propagation, solidification and multiphase flows. Another problem involving moving boundaries is that of fluid-structure interactions, and I am currently also working on simulating a suspension of flexible fibers immersed in a Newtonian fluid.

One type of numerical methods for moving interface problems are the so called interface tracking methods, that are developed to numerically describe and track such interfaces when they move, deform and interact as determined by the underlying physics. Initially, I worked with a front-tracking method and later a level-set method, coupled to a finite element discretization of the Navier-Stokes equation, to simulate multiphase flows. In my PhD thesis, I also developed a new interface tracking technique, the segment projection method. This method can be viewed as a compromise between the front-tracking and level-set methods. It relies on the partitioning of an interface into several parts. In the case of two dimensional calculations, the interfaces are curves in the plane. Each curve is described by a union of overlapping curve segments, where the segments are chosen such that they can be given as functions of one spatial variable. The segment functions are discretized using one-dimensional Eulerian grids in the x and y directions. Hence, the curve is explicitly discretized with points on the curve as in the front-tracking method, thereby keeping the lower dimensionality, but the discretization is Eulerian as in the level-set method. The Eulerian form makes it easier to handle topological changes as in the example of merging bubbles.

The first application of the segment projection method was the multiphase flow simulations that I had earlier applied both the front-tracking and the level-set methods to. In this context, also the segment projection method is integrated into a finite element approximation of the incompressible two-dimensional Navier-Stokes equations. The segment projection method has also later been applied to front propagation in geometrical optics and to simulations of epitaxial growth (applicable to growth of thin films in semi conductor production). This work is done in collaboration with Björn Engquist, and for the geometrical optics application also with Olof Runborg. In both these applications, the Eulerian explicit discretization of the interfaces is exploited, as we solve partial differential equations on the moving front/interface.

Very frequently in moving interface problems, there are forces acting on the interfaces, such as surface tension forces, elastic forces etc. When using an interface tracking method to couple the interface dynamics to for example the fluid equations to be solved on the background grid, these forces with singular support must be discretized on that grid. Motivated by this, I have been studying how singular forces and discontinuous coefficients should be discretized. First, in connection to quadrature, as motivated by a finite element discretizations, and later also regarding finite difference discretizations of different differential equations.

Dynamics of multiple interacting flexible fibers in Stokes flow

We have developed a numerical method for simulating slender flexible fibers suspended in fluids. This is work done in collaboration with M. Shelley.

Here are some snapshots from a simulation of 25 flexible fibers in a shear flow. The fibers are colored by the tension in the fibers.
<

The understanding of the dynamics of such suspensions are fundamental to understanding many flows arising in physics, biology and engineering. Examples include fiber-reinforced composites, the dynamics and rheology of biological polymers and the motility of microscopic organisms. Such filaments often have aspect ratios of length to radius ranging from a few hundred to several thousand. Full discretizations of such thin objects in a 3D domain is very costly. The Reynolds numbers for the applications listed are however very low.

We consider a three-dimensional flow with multiple, interacting fibers and make explicit use of both the fact that we are considering Stokes flow, as well as of the slenderness of the fibers. The key points are that for Stokes flow, boundary integral methods can be employed to reduce the three-dimensional dynamics to the dynamics of the two-dimensional filament surfaces, and that using slender body asymptotics, this can be further reduced to the dynamics of the one-dimensional filament centerlines. The resulting intergral equations take into account the fluid-filament interaction as well as filament-filament interactions, as mediated by the fluid.

The filaments we are considering are inextensible and elastic. Replacing the force in the integral equation by an explicit expression based on the shape of the filament, using Euler-Bernoulli elasticity, yields a time-dependent integral equation for the motion of the filament centerline. This equation will include the filament tension, for which an equation can be derived by using the condition of inextensibility of the filament. The evolution of the system therefore, in each time step, requires the solution of an auxiliary integro-differential equation for the filament tension.

In our numerical algorithm, the filament centerlines are parameterized by arclength, and discretized uniformly. All derivatives are computed using second-order divided differences. Special quadrature rules have been developed to compute the necessary integrals. A second-order time-stepping scheme is used, with implicit treatment of high derivatives. A product integration method is applied to evaluate the integrals in the integral operators.

We place our filaments in a shear flow. In the non-dimensionalized equations, there are two physical parameters. One is the ratio of the radius to the length, and the other one is a parameter Z that relates the characteristic fluid drag to the filament elasticity. If one single straight filament is inserted in a plane shear flow at some angle to the x-axis, it will simply rotate around its center of mass. If we introduce a small perturbation to the filament, so that it is not exactly straight, there are two possible scenarios. Either, this perturbation will disappear with time, and the filament will stay straight. However, if the filament is under compression for some time, and if the value of Z is large enough, then a so called buckling occurs. As we increase Z, this buckling becomes more pronounced.

Next follows three short animations that show the fundamentals for this buckling, for a length to radius ration of 1000, for three different values of Z. In these animations, the filament has been colored with the line tension. The line tension changes from negative to positive, as the filament goes from being under compression to being under extension.

Each movie is approx. 700KBytes

Next animation is one of 25 interacting filaments, inserted into a oscillating background shear flow. We impose periodic boundary conditions in the streamwise (x) direction. The animation shows one period in time. (Filaments not colored by line-tension in this movie).

Movie is approx. 2Mbytes

Rigid fiber suspensions

Up to A-K. Tornberg's home page.