I hope that all of my readers have had a restful Thanksgiving Break and are ready to begin their final exam studies. This week, we will discuss how to design a numerical program that will predict the path a particle will take in a real-world fluid flow. To do this, we must first outline our overall problem.
Suppose that we have mounted a wing inside a wind tunnel. When we turn the tunnel on, air will blow past the wing. This special air will have tiny particles distributed within it that will travel along streamlines as the air moves past the wing. We will then use a camera to take a picture of the flow around the wing after time intervals of 50 milliseconds. Because the tiny particles are present in the air, we will be able to tell from the photographs how the air moves past the wing. Using a computer program in Matlab, we will use the photos to generate a discretized velocity vector field, U(x,y,z,t) for the air moving around the wing. In this situation, U has three components, each of which depends on our three spatial variables and time.
U(x,y,z,t) = <u(x,y,z,t), v(x,y,z,t), w(x,y,z,t)>
Now suppose that in an upcoming experiment, we want to track the motion of one water molecule as it moves through the 3-dimensional flow described by U(x,y,z,t). This molecule will start out with spatial coordinates of (x0, y0, z0) at time t = t0 and before our next experiment, we want to predict the molecule’s spatial coordinates at any future time t.
For the sake of compactness, let us write X0 = <x0, y0, z0> and X(t) = <x(t), y(t), z(t)>, where X(t) is a vector-valued function for the position of our molecule and X0 is a vector representing the molecule’s initial position.
The Second Fundamental Theorem of Calculus applied to vector-valued functions gives us the following for any time t after t0. (The integral runs from t0 to t.)
X(t) = X0 + ∫ (dX/dt) dt
When taking the integral of a vector-valued function, we integrate each of the function’s components separately. It also might seem a little confusing to be integrating a function of t from t0 to t. What we are doing here is finding an antiderivative for dX/dt and then evaluating that antiderivative at t = t, and then subtracting away the value of the antiderivative at t = t0.
Furthermore, we know that the velocity of the molecule, dX/dt, is equal to U(x,y,z,t), giving us a system of three coupled differential equations.
dx/dt = u(x,y,z,t)
dy/dt = v(x,y,z,t)
dz/dt = w(x,y,z,t)
Compactly, we can write this entire system as dX/dt = U(x,y,z,t), thus yielding:
X(t) = X0 + ∫ U(x,y,z,t) dt, where the integral is from time t0 to time t.
The complexity in solving this problem comes in evaluating the integral above. To do this accurately, we would need to guess what U(x,y,z,t) would be at positions and times for which our photo analysis program has not already given us a value of U. Moreover, we must choose an algorithm for numerically solving this integral that is both fast and accurate. In my own work, I have adapted what is called a Runge-Kutta-Fehlberg 45 (RKF 45) integration scheme to predict the position of a molecule moving through U. While the details of this process are a bit too detailed for one blog post, I urge readers to learn more about this method through Wikipedia or YoutTube videos. Determining particle trajectories can be particularly useful because they have allowed engineers to then figure out how quickly two particle paths can diverge away from each other. Once they have this information, they could gain a much better understanding of the overall structure of the flow they are dealing with (the interested reader can study Finite-Time Lyapunov Exponents and Lagrangian Coherent Structures to learn more).
In the next blog post, we will turn our attention to the theory of electrostatics and Poisson’s equation (∂2u/∂x2 + ∂2u/∂y2 = f(x,y)). Until then, take care.