1
\documentclass[12pt]{article}
4
\usepackage{amsmath,amssymb}
8
\title{Two Dimensional Euler Simulation}
11
S.J. Montgomery-Smith\\
12
Department of Mathematics\\
13
University of Missouri\\
14
Columbia, MO 65211, U.S.A.\\
15
stephen@math.missouri.edu\\
16
http://www.math.missouri.edu/\~{}stephen}
18
\date{September 10, 2000}
22
This document describes a program I wrote to simulate the
23
two dimensional Euler Equation --- a program that is part
24
of the {\tt xlock} screensaver as the {\tt euler2d}
25
mode. A similar explanation may also be found in the
26
book by Chorin \cite{C}.
28
\section{The Euler Equation}
30
The Euler Equation describes the motion of an incompressible
31
fluid that has no viscosity. If the fluid is contained
32
in a domain $\Omega$ with boundary $\partial \Omega$, then
33
the equation is in the vector field $u$ (the velocity)
35
scalar field $p$ (the pressure):
37
\frac{\partial}{\partial t} u &=& -u \cdot \nabla u + \nabla p \\
38
\nabla \cdot u &=& 0 \\
39
u \cdot n &=& 0 \quad \text{on $\partial \Omega$}
41
where $n$ is the unit normal to $\partial \Omega$.
45
It turns out that it can be easier write these equations
46
in terms of the vorticity. In two dimensions the vorticity
47
is the scalar $w = \partial u_2/\partial x - \partial u_1/\partial y$.
48
The equation for vorticity becomes
49
\[ \frac{\partial}{\partial t} w = -u \cdot \nabla w .\]
50
A solution to this equation can be written as follows. The velocity
51
$u$ causes a flow, that is, a function $\varphi(t,x)$ that tells where
52
the particle initially at $x$ ends up at time $t$, that is
54
\frac\partial{\partial t} \varphi(t,x)
55
= u(t,\varphi(t,x)) .\]
57
for $w$ tells us that the vorticity is ``pushed'' around by the flow,
58
that is, $w(t,\varphi(t,x)) = w(0,x)$.
60
\section{The Biot-Savart Kernel}
62
Now, once we have the vorticity, we can recover the velocity $u$ by
65
\partial u_2/\partial x - \partial u_1/\partial y &=& w \\
66
\nabla \cdot u &=& 0 \\
67
u \cdot n &=& 0 \quad \text{on $\partial \Omega$}.
69
This equation is solved by using a Biot-Savart kernel $K(x,y)$:
70
$$ u(x) = \int_\Omega K(x,y) w(y) \, dy .$$
71
The function $K$ depends upon the choice of domain. First let us consider
72
the case when $\Omega$ is the whole plane (in which case the boundary
73
condition $u \cdot n = 0$ is replaced by saying that $u$ decays at infinity).
76
K(x,y) = K_1(x,y) = c \frac{(x-y)^\perp}{|x-y|^2} .
78
Here $x^\perp = (-x_2,x_1)$, and $c$ is a constant, probably something
79
like $1/2\pi$. In any case we will set it to be one, which in effect
80
is rescaling the time variable, so we don't need to worry about it.
82
We can use this as a basis to find $K$ on the unit disk
83
$\Omega = \Delta = \{x:|x|<1\}$. It turns out to be
85
K_2(x,y) = K_1(x,y) - K_1(x,y^*) ,
87
where $y^* = y/|y|^2$ is called the reflection of $y$ about the
88
boundary of the unit disk.
90
Another example is if we have a bijective analytic function
91
$p:\Delta \to {\mathbb C}$, and we let $\Omega = p(\Delta)$.
92
(Here we think of $\Delta$ as a subset of $\mathbb C$, that is,
93
we are identifying the plane with the set of complex numbers.)
95
\[ K_p(p(x),p(y)) = K_2(x,y)/|p'(x)|^2 .\]
96
Our simulation considers the last case. Examples of such
97
analytic functions include series
98
$p(x) = x + \sum_{n=2}^\infty c_n x^n$, where
99
$\sum_{n=2}^\infty n |c_n| \le 1$.
100
(Thanks to David Ullrich for pointing this out to me.)
102
\section{The Simulation}
104
Now let's get to decribing the simulation. We assume a rather
105
unusual initial distribution for the vorticity --- that the
106
vorticity is a finite sum of dirac delta masses.
107
\[ w(0,x) = \sum_{k=1}^N w_k \delta(x-x_k(0)) .\]
108
Here $x_k(0)$ is the initial place where the points
109
of vorticity are concentrated, with values $w_k$.
110
Then at time $t$, the vorticity becomes
111
\[ w(t,x) = \sum_{k=1}^N w_k \delta(x-x_k(t)) .\]
112
The points of fluid $x_k(t)$ are pushed by the
113
flow, that is, $x_k(t) = \varphi(t,x_k(0))$, or
114
\[ \frac{\partial}{\partial t} x_k(t) = u(t,x_k(t)) .\]
115
Putting this all together, we finally obtain the equations
116
\[ \frac{\partial}{\partial t} x_k = \alpha_k \]
118
\[ \alpha_k = \sum_{l=1}^N w_l K(x_k,x_l) .\]
119
This is the equation that our simulation solves.
121
In fact, in our case, where the domain is $p(\Delta)$,
122
the points are described by points
123
$\tilde x_k$, where $x_k = p(\tilde x_k)$. Then
127
\frac{\partial}{\partial t} \tilde x_k &=& \tilde\alpha_k \\
129
\tilde\alpha_k &=& \frac1{|p'(\tilde x_k)|^2}
130
\sum_{l=1}^N w_l K_2(\tilde x_k,\tilde x_l) .
133
We solve this $2N$ system of equations using standard
134
numerical methods, in our case, using the second order midpoint method
135
for the first step, and thereafter using the second order Adams-Bashforth
136
method. (See for example the book
137
by Burden and Faires \cite{BF}).
139
\section{The Program - Data Structures}
141
The computer program solves equation (\ref{tildex-p1}), and displays
142
the results on the screen, with a boundary. All the information
143
for solving the equation and displaying the output is countained
144
in the structure {\tt euler2dstruct}. Let us describe some of
145
the fields in {\tt euler2dstruct}.
146
The points $\tilde x_k$ are contained
147
in {\tt double *x}: with the coordinates of
148
$\tilde x_k$ being the two numbers
149
{\tt x[2*k+0]}, {\tt x[2*k+1]}. The values $w_k$ are contained
150
in {\tt double *w}. The total number of points is
151
{\tt int N}. (But only the first {\tt int Nvortex} points
152
have $w_k \ne 0$.) The coefficients of the analytic function
153
(in our case a polynomial) $p$
154
are contained in {\tt double p\_coef[2*(deg\_p-1)]} --- here
155
{\tt deg\_p} is the degree of $p$, and the real and imaginary
156
parts of the coefficient
157
$c_n$ is contained in {\tt p\_coef[2*(n-2)+0]} and {\tt p\_coef[2*(n-2)+1]}.
159
\section{Data Initialization}
161
The program starts in the function {\tt init\_euler2d}. After allocating
162
the memory for the data, and initialising some of the temporary variables
163
required for the numerical solving program, it randomly assigns the
164
coefficients of $p$, making sure that $\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$.
165
Then the program figures out how to draw the boundary, and what rescaling
166
of the data is required to draw it on the screen. (This uses the
167
function {\tt calc\_p} which calculates $p(x)$.)
169
Next, it randomly assigns the initial values of $\tilde x_k$. We want
170
to do this in such a way so that the points are uniformly spread over the
171
domain. Let us first consider the case when the domain is the unit circle
172
$\Delta$. In that case the proportion of points that we would expect
173
inside the circle of radius $r$ would be proportional to $r^2$. So
175
\[ r = \sqrt{R_{0,1}},\quad \theta = R_{-\pi,\pi}, \quad
176
\tilde x_k = r (\cos \theta, \sin \theta) .\]
177
Here, and in the rest of this discussion, $R_{a,b}$ is a function
178
that returns a random variable uniformly distributed over the interval
181
This works fine for $\Delta$, but for $p(\Delta)$, the points
182
$p(\tilde x_k)$ are not uniformly distributed over $p(\Delta)$,
183
but are distributed with a density proportional to
184
$1/|p'(\tilde x_k)|^2$. So to restore the uniform density we need
185
to reject this value of $\tilde x_k$ with probability proportional
186
to $|p'(\tilde x_k)|^2$. Noticing that the condition
187
$\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$ implies that
188
$|p'(\tilde x_k)| \le 2$, we
189
do this by rejecting if $|p'(\tilde x_k)|^2 < R_{0,4}$.
190
(This makes use of the function {\tt calc\_mod\_dp2} which calculates
193
\section{Solving the Equation}
195
The main loop of the program is in the function {\tt draw\_euler2d}.
196
Most of the drawing operations are contained in this function, and
197
the numerical aspects are sent to the function {\tt ode\_solve}.
198
But there is an aspect of this that I would like
199
to discuss in the next section, and so we will look at a simple method for
200
numerically solving differential equations.
203
(nothing to do with the Euler Equation), is as
204
follows. Pick a small number $h$ --- the time step (in
205
the program call {\tt delta\_t}). Then we approximate
206
the solution of the equation:
208
\label{method-simple}
209
\tilde x_k(t+h) = \tilde x_k(t) + h \tilde\alpha_k(t) .
211
The more sophisticated methods we use are variations of
212
the Euler Method, and so the discussion in the following section
215
In the program, the quantities $\tilde\alpha_k$, given by
216
equations (\ref{tildex-p2}) are calculated by the function
218
(which in turns calls {\tt calc\_all\_mod\_dp2} to
219
calculate $|p'(\tilde x_k)|^2$ at all the points).
222
\section{Subtle Perturbation}
224
Added later: the scheme described here seems to not be that effective,
225
so now it is not used.
227
One problem using a numerical scheme such as the Euler Method occurs
228
when the points $\tilde x_k$ get close to the boundary
229
of $\Delta$. In that case, it is possible that the new
230
points will be pushed outside of the boundary. Even if they
231
are not pushed out of the boundary, they may be much closer
232
or farther from the boundary than they should be.
233
Our system of equations is very sensitive to how close points
234
are to the boundary --- points with non-zero vorticity
235
(``vortex points'') that are close to the boundary travel
236
at great speed alongside the boundary, with speed that is
237
inversely proportional to the distance from the boundary.
239
A way to try to mitigate this problem is something that I call
240
``subtle perturbation.''
242
the unit disk to points in the plane using the map
244
F(x) = f(|x|) \frac x{|x|} ,
246
where $f:[0,1]\to[0,\infty]$ is an increasing continuous
247
bijection. It turns out that a good choice is
251
(The reason for this is that points close to each other
253
about $r$ from the boundary will be pushed around so that
254
their distance from each other is about multiplied by the
255
derivative of $\log r$, that is, $1/r$.)
256
Note that the inverse of this function is given by
258
F^{-1}(x) = f^{-1}(|x|) \frac x{|x|} ,
262
f^{-1}(t) = 1-e^{-t} .
265
So what we could do is the following: instead of working with
266
the points $\tilde x_k$, we could work instead with the points
267
$y_k = F(\tilde x_k)$. In effect this is what we do.
268
Instead of performing the computation (\ref{method-simple}),
269
we do the calculation
271
y_k = F(\tilde x_k(t)) + h {\cal A}(\tilde x_k) \tilde\alpha_k(t)
274
${\cal A}(x)$ is the matrix of partial derivatives of $F$:
285
\left(\frac{f'(|x|)}{|x|} - \frac{f(|x|)}{|x|^2}\right)
288
x_{1}^2 & x_{1} x_{2}\\
289
x_{1} x_{2} & x_{2}^2
295
\tilde x_k(t+h) = F^{-1}(y_k).
297
These calculations are done in the function {\tt perturb}, if
298
the quantity {\tt SUBTLE\_PERTURB} is set.
300
\section{Drawing the Points}
302
As we stated earlier, most of the drawing functions are contained
303
in the function {\tt draw\_euler2d}. If the variable
304
{\tt hide\_vortex} is set (and the function {\tt init\_euler2d}
305
will set this with probability $3/4$), then we only display
306
the points $\tilde x_k$ for ${\tt Nvortex} < k \le N$. If
307
{\tt hide\_vortex} is not set, then the ``vortex points''
308
$\tilde x_k$ ($1 \le k \le {\tt Nvortex}$) are displayed in white.
309
In fact the points $p(\tilde x_k)$ are what are put onto the screen,
310
and for this we make use of the function {\tt calc\_all\_p}.
312
\section{Addition to Program: Changing the Power Law}
314
A later addition to the program adds an option {\tt eulerpower},
315
which allows one to change the power law that describes how
316
the vortex points influence other points. In effect, if this
317
option is set with the value $m$, then the Biot-Savart Kernel
319
$$ K_1(x,y) = \frac{(x-y)^\perp}{|x-y|^{m+1}}, $$
321
$$ K_2(x,y) = K_1(x,y) - |y|^{1-m} K_1(x,y) .$$
322
So for example, setting $m=2$ corresponds to the
323
quasi-geostrophic equation. (I haven't yet figured out
324
what $K_p$ should be, so if $m \ne 1$ we use the unit circle
327
\begin{thebibliography}{9}
329
\bibitem{BF} Richard L. Burden, J. Douglas Faires, Numerical Analysis,
330
sixth edition, Brooks/Cole, 1996.
332
\bibitem{C} Alexandre J. Chorin, Vorticity and Turbulence,
333
Applied Mathematical Sciences, Vol 103, Springer Verlag, 1994.
335
\end{thebibliography}
added
removed
hacks/euler2d.tex
