Flow Simulation

Here we present the simulation results of flows around a number of rigid and variable shape bodies. The different cases of Potential, Stokesian, and Navier-Stokes flow are considered.

The linear Potential and Stokes flows are solved using a simple panel method, whereas the Navier-Stokes flow is solved using the Immersed Boundary Method.

Pulsating Circle

We start with the simple example of a pulsating circle. It’s boundary is described by:

\[r(\theta) = 1+q_0(t)\]

with the periodic shape function

\[q_0(t) = \sin(t)\]

In python code this is expressed as

shapevariables = (
        lambda t: sin(t),

X = lambda s,q: cos(s) (1+q[0])
Y = lambda s,q: sin(s) (1+q[0])
bodies = (

The solution can easily be found from the continuity equaiton alone

\[u_r(\theta) \propto \frac{1}{r}\]

and the numerical simulation is run just for validation purposes. Whereas the panel methods work fine in this case, the immersed boundary method does not converge as it has not been designed to handle cases where the panel size on the boundary changes.

Below is the result of both the stokes and potential solver (the result is the same for both)


Translating Circle

We now consider the periodic translation of a rigid circle. Its boundary is described by:

\[\begin{split}x(\theta) &= \cos(\theta) + q_0(t)\\ y(\theta) &= \sin(\theta)\end{split}\]

with a periodic shape function

\[q_0(t) = \frac{1}{2}\sin(t)\]

The code listing

shapevariables = (
        lambda t: 0.5*sin(t),

X = lambda s,q: cos(s) + q[0]
Y = lambda s,q: sin(s)
bodies = (

Due to Stoke’s Paradox there does not exist a 2D creeping flow around a translating cylinder that vanishes at infinity. The Stokesian solver handles this problem gracefully and outputs a uniform flow for \(Re=0\). The immersed boundary method is used to solve for the \(Re=20\), \(Re=100\) and \(Re=1000\) flows. And the potential solver (assumes irrotational flow) is used to solve for \(Re=\infty\). The results are presented in the table below:

Re=0 Re=1000
im1 im4
Re=20 Potential, irrotatoinal
im2 im5

There is an inverse correlation between boundary thickness and Reynolds number as expected

\[\delta \propto \frac{1}{\sqrt{Re}}\]


You can view the simulations at higher resolution by Right Click->View Image (on Firefox) or using an equivalent method on another browser.

Mason & Burdick Amoeba

Here we take the example given by [MB99]. The boundary is simulated by:

\[r(\theta) = 1 + \epsilon (s_1(t) \cos(2\theta) + s_2(t) \sin(3\theta) )\]

By using an appropriate time evolution of the two shape modes we can achieve forward motion.


In the simulations below, the parameter is set to \(\epsilon=0.1\) and the following evolution in time of the shape modes is used

\begin{eqnarray*} s_{1}(t)&=&1-\cos t\\ s_{2}(t)&=&\sin t \end{eqnarray*}

In order for the immersed boundary method to converge we must also ensure that the area enclosed by the boundary remains constant. This can be achieved by adding a normalization:

\[r(\theta) = \frac{1}{N(t)}\big(1 + \epsilon (s_1(t) \cos(2\theta) + s_2(t) \sin(3\theta) )\big)\]


\[N(t) = \sqrt {1 + \frac{1}{2}\epsilon^2\left ( {s_1(t)}^2 + {s_2(t)}^2 \right)}\]


shapevariables = (
        lambda t: 1-cos(t),
        lambda t: sin(t),

epsilon = 0.1
R = lambda s,q: 1 + epsilon*(q[0]*cos(2*s) + q[1]*cos(3*s))
N = lambda s,q: sqrt(1+0.5*epsilon**2*(q[0]**2+q[1]**2))
X = lambda s,q: R(s,q)*cos(s)/N(s,q)
Y = lambda s,q: R(s,q)*sin(s)/N(s,q)
bodies = (
Re=0 Re=1000
Re=10 Potential, irrotatoinal

Shapere & Wilczek Amoeba

Here we take a different exmample of an amoeba studied by Shapere and Wilczek in [SW89].

This time the boundary is given by:

\begin{eqnarray*} x(\theta)&=&\cos\theta\\ &&\quad+s_{1}(t)(a\cos\theta+b\sin2\theta)\\ &&\quad+s_{2}(t)(-a\cos2\theta+b\sin\theta)\\ y(\theta)&=&\sin\theta\\ &&\quad+s_{1}(t)(-a\sin\theta+b\cos2\theta)\\ &&\quad+s_{2}(t)(a\sin2\theta+b\cos\theta) \end{eqnarray*}

Below we use the shape modes evolving in time as follows

\begin{eqnarray*} s_{1}(t)&=&\cos t\\ s_{2}(t)&=&\sin t \end{eqnarray*}

and the parameters are set to \(a=\frac{3}{10}\quad b=\frac{3}{200}\)

We will not compute the immersed boundary case here, as ensuring the area remains constant in this case is not trivial.

Re=0 Potential, irrotation


[MB99]Richard Mason and Joel Burdick. Propulsion and control of deformable bodies in an ideal fluid. In Robotics and Automation, 1999. Proceedings. 1999 IEEE International Conference on, volume 1, 773–780. IEEE, 1999.
[SW89]Alfred Shapere and Frank Wilczek. Geometry of self-propulsion at low reynolds number. Journal of Fluid Mechanics, 198:557–585, 1989.