Lateral advection

PyLag can be easily used with analytical models of different flow fields. This provides a convenient means of testing the model implementation. The example described here is taken from Kreyszig, E. (2006) Advanced Engineering Mathematics, Ch. 18. P762.

Background theory

In the example, the complex potential with \(F(z) = z^{2} = x^{2} - y^{2} + 2ixy\) models a flow with equipotential lines and streamlines given by:

\[\begin{split}\begin{eqnarray} \Phi &=& x^2 - y^2 &=& const \\ \Psi &=& 2xy &=& const \end{eqnarray}\end{split}\]

The velocity vector is:

\[\begin{equation} \mathbf{u} = 2(x - iy) \end{equation}\]

with components:

\[\begin{split}\begin{eqnarray} u &=& &2&x \\ v &=& -&2&y \end{eqnarray}\end{split}\]

The speed of the flow is given by:

\[\begin{equation} |\mathbf{u}| = \sqrt{x^2 + y^2} \end{equation}\]

The flow may be interpreted as the flow in a channel bounded by the positive coordinate axes and a hyperbola, given by \(xy = const\). To calculate particle trajectiories, PyLag solves the equation:

\[\begin{equation} \label{eqn:partmotion} \dfrac{\partial}{\partial t} \mathbf{X}_{i}(t, \mathbf{r}_{i}) = \mathbf{U}_{i}(t, \mathbf{X}_{i}) \end{equation}\]

where \(\mathbf{r}_{i} = \mathbf{X}_i(t=t_0)\) is the position vector of particle \(i\) at time \(t=t_{0}\), and \(\mathbf{U}_{i}\) is the particle’s velocity vector. Here, \(\mathbf{U}_{i} = [\mathbf{u}(t, \mathbf{x})]_{\mathbf{x} = \mathbf{X}_{i}}\) where \(\mathbf{u}\) is the fluid velocity vector, as given above.

In practise, PyLag computes solutions to the above equation using a Random Displacement Model of the form:

\[\begin{equation} dX_{j} = \left[u_{j}\left(\mathbf{x},t\right) + \dfrac{\partial D_{jk}\left(\mathbf{x},t\right)}{\partial x_{k}}\right]dt + \left(2 D_{jk}\left(\mathbf{x},t\right)\right)^{1/2}dW_{k}(t) \end{equation}\]

Here, \(dX_{j} = d\mathbf{X}\) is the incremental change in a particle’s position and \(D_{jk}\) is the diffusion coefficient or tensor. The first term on the right-hand side (RHS) of the equation is a deterministic drift term; the second is a stochastic term that models the action of diffusion. \(dW_{j}(t)\) is an incremental Wiener process that builds stochasticity into the model.

In this particular problem, the diffusion coefficient is zero. In component form, the model then reduces to:

\[\begin{equation} \label{eqn:dx} dX = u\left(\mathbf{x},t\right) dt, \end{equation}\]

\[\begin{equation} \label{eqn:dy} dY = v\left(\mathbf{x},t\right) dt, \end{equation}\]

\[\begin{equation} \label{eqn:dz} dZ = 0, \end{equation}\]

which can be solved using one of the standard numerical integration schemes supported by PyLAg.

Running the particle tracking model

To compute particle trajectories within the given flow field using PyLag, we first subclass PyLag’s DataReader. The new subclass must be implemented in Cython, since several of its functions have arguments without a direct counterpart in Python. The above model has been implemented in pylag.mock, and is called MockVelocityDataReader.

Note:

It is possible to extend PyLag by implementing alternative concrete DataReader subclasses in pylag.mock. Alternatively, you may choose to implement these in a separate module. Either way, each module must first by cythonized and compiled before it can be used.

We will start three particles off with initial coordinates \(X_1(t = t_0) = (0.125, 1.)\), \(X_2(t = t_0) = (0.3525, 1.)\) and \(X_3(t = t_0) = (0.5, 1.)\). The corresponding streamlines are given by \(\Psi = \{0.5, 1, 2\}\). Given the initial particle coordinates, we can use PyLag to compute particle pathlines, which can be compared directly with the above streamlines. To do this, we use one of PyLag’s numerical integration schemes. The integration is managed by an object of type pylag.numerics.NumIntegrator, which in turn farms out most of the hard work to an object of type pylag.numerics.ItMethod. Given the type of flow field, we will use StdNumIntegrator to manage the integration, and the Adv_RK4_2D iterative method, which encodes a fourth order runga kutta scheme. The following code uses PyLag to compute particle positions as a function of time:

[1]:

import warnings
import numpy as np
import matplotlib
from configparser import SafeConfigParser

from pylag.numerics import StdNumMethod
from pylag.particle_cpp_wrapper import ParticleSmartPtr
from pylag.mock import MockVelocityDataReader

# Ensure inline plotting
%matplotlib inline

warnings.filterwarnings('ignore')

# Set the time step to 0.1s
time_step = 0.01

# PyLag reads many of its parameters from a run config, which we create here
config = SafeConfigParser()
config.add_section('SIMULATION')
config.add_section('NUMERICS')

# Set the coordinate system, which here is a simple cartesian coordinate system
config.set('SIMULATION', 'coordinate_system', 'cartesian')

# Set the type of iterative method to a 2D fourth order Runga Kutta method
config.set('NUMERICS', 'iterative_method', 'Adv_RK4_2D')

# Set the time step
config.set('NUMERICS', 'time_step_adv', str(time_step))

# The data reader
data_reader = MockVelocityDataReader()

# The numerical integration scheme - this uses the config to initialise its state
num_method = StdNumMethod(config)

# The three particles, saved in a dictionary.
particles = {'Particle1': ParticleSmartPtr(x1=0.125, x2=1.0),
             'Particle2': ParticleSmartPtr(x1=0.25, x2=1.0),
             'Particle3': ParticleSmartPtr(x1=0.375, x2=1.0)}

# Dictionary in which particle position data will be saved. x and y positions are stored
# in lists. Initialise this with the particle's position at t = t_0.
positions_out = {}
for key in particles.keys():
    positions_out[key] = {'x1': [particles[key].x1], 'x2': [particles[key].x2]}

# We will start the model off at t = 0 s and integrate it forward for 2 s
time_grid = np.arange(0.0, 2.0, time_step).tolist()

# Perform the integration
for time in time_grid:
    for key in particles.keys():
        # Udate partcle positions.
        num_method.step_wrapper(data_reader, time, particles[key])

        # Save particle positions
        positions_out[key]['x1'].append(particles[key].x1)
        positions_out[key]['x2'].append(particles[key].x2)

Visualising the result

Below we animate the particle positions on top of the streamlines. A quiver plot of the velocity field is also shown.

[2]:

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

from matplotlib import animation, rc
from IPython.display import HTML

font_size=15

fig, ax = plt.subplots(figsize = (8, 8))

# Plot the line x = y
ax.plot([0,1], [0,1], 'k--')

# Add a quiver plot of the velocity field
x_grid, y_grid = np.meshgrid(np.linspace(0,1,20),np.linspace(0,1,20))
quiver_plot = ax.quiver(x_grid, y_grid, 2.*x_grid, -2.*y_grid)

# Add streamlines
y = np.linspace(0.01, 5., 100, dtype=float)
for psi in [0.25, 0.5, 0.75]:
    x = 0.5 * psi / y

    ax.plot(x, y, c='b', zorder=1)

# Configure the plot
ax.set_xlim((0, 1))
ax.set_ylim((0, 1))

ax.set_xlabel('x', fontsize=font_size)
ax.set_ylabel('y', fontsize=font_size)

ax.tick_params(axis='both', which='major', labelsize=font_size)
ax.set_aspect('equal')

# Animate particle trajectories
scatter = ax.scatter([], [], c='r', s=100, zorder=2)

def animate(i):
    xpos = []
    ypos = []
    for key in particles.keys():
        xpos.append(positions_out[key]['x1'][i])
        ypos.append(positions_out[key]['x2'][i])

    data = np.array([xpos, ypos])
    scatter.set_offsets(data.transpose())
    return (scatter,)

# Prevent the basic figure from being plotted too
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, frames=110, interval=20, blit=True)

HTML(anim.to_html5_video())

[2]:

It is clear that all three particles closely follow their respective streamlines, as they should. The dashed line marks the point along each streamline where the speed is at a minimum. This can be seen in the particle trajectories, with particles decelerating as they approach the point where the dashed line intersects each streamline, then accelerating away from it.