338 lines
12 KiB
TeX
338 lines
12 KiB
TeX
|
\documentclass[12pt]{article}
|
||
|
|
|
||
|
|
%\usepackage{fullpage}
|
||
|
|
\usepackage{amsmath,amssymb}
|
||
|
|
|
||
|
|
\begin{document}
|
||
|
|
|
||
|
|
\title{Two Dimensional Euler Simulation}
|
||
|
|
|
||
|
|
\author{
|
||
|
|
S.J. Montgomery-Smith\\
|
||
|
|
Department of Mathematics\\
|
||
|
|
University of Missouri\\
|
||
|
|
Columbia, MO 65211, U.S.A.\\
|
||
|
|
stephen@math.missouri.edu\\
|
||
|
|
http://www.math.missouri.edu/\~{}stephen}
|
||
|
|
|
||
|
|
\date{September 10, 2000}
|
||
|
|
|
||
|
|
\maketitle
|
||
|
|
|
||
|
|
This document describes a program I wrote to simulate the
|
||
|
|
two dimensional Euler Equation --- a program that is part
|
||
|
|
of the {\tt xlock} screensaver as the {\tt euler2d}
|
||
|
|
mode. A similar explanation may also be found in the
|
||
|
|
book by Chorin \cite{C}.
|
||
|
|
|
||
|
|
\section{The Euler Equation}
|
||
|
|
|
||
|
|
The Euler Equation describes the motion of an incompressible
|
||
|
|
fluid that has no viscosity. If the fluid is contained
|
||
|
|
in a domain $\Omega$ with boundary $\partial \Omega$, then
|
||
|
|
the equation is in the vector field $u$ (the velocity)
|
||
|
|
and the
|
||
|
|
scalar field $p$ (the pressure):
|
||
|
|
\begin{eqnarray*}
|
||
|
|
\frac{\partial}{\partial t} u &=& -u \cdot \nabla u + \nabla p \\
|
||
|
|
\nabla \cdot u &=& 0 \\
|
||
|
|
u \cdot n &=& 0 \quad \text{on $\partial \Omega$}
|
||
|
|
\end{eqnarray*}
|
||
|
|
where $n$ is the unit normal to $\partial \Omega$.
|
||
|
|
|
||
|
|
\section{Vorticity}
|
||
|
|
|
||
|
|
It turns out that it can be easier write these equations
|
||
|
|
in terms of the vorticity. In two dimensions the vorticity
|
||
|
|
is the scalar $w = \partial u_2/\partial x - \partial u_1/\partial y$.
|
||
|
|
The equation for vorticity becomes
|
||
|
|
\[ \frac{\partial}{\partial t} w = -u \cdot \nabla w .\]
|
||
|
|
A solution to this equation can be written as follows. The velocity
|
||
|
|
$u$ causes a flow, that is, a function $\varphi(t,x)$ that tells where
|
||
|
|
the particle initially at $x$ ends up at time $t$, that is
|
||
|
|
\[
|
||
|
|
\frac\partial{\partial t} \varphi(t,x)
|
||
|
|
= u(t,\varphi(t,x)) .\]
|
||
|
|
Then the equation
|
||
|
|
for $w$ tells us that the vorticity is ``pushed'' around by the flow,
|
||
|
|
that is, $w(t,\varphi(t,x)) = w(0,x)$.
|
||
|
|
|
||
|
|
\section{The Biot-Savart Kernel}
|
||
|
|
|
||
|
|
Now, once we have the vorticity, we can recover the velocity $u$ by
|
||
|
|
solving the equation
|
||
|
|
\begin{eqnarray*}
|
||
|
|
\partial u_2/\partial x - \partial u_1/\partial y &=& w \\
|
||
|
|
\nabla \cdot u &=& 0 \\
|
||
|
|
u \cdot n &=& 0 \quad \text{on $\partial \Omega$}.
|
||
|
|
\end{eqnarray*}
|
||
|
|
This equation is solved by using a Biot-Savart kernel $K(x,y)$:
|
||
|
|
$$ u(x) = \int_\Omega K(x,y) w(y) \, dy .$$
|
||
|
|
The function $K$ depends upon the choice of domain. First let us consider
|
||
|
|
the case when $\Omega$ is the whole plane (in which case the boundary
|
||
|
|
condition $u \cdot n = 0$ is replaced by saying that $u$ decays at infinity).
|
||
|
|
Then
|
||
|
|
\begin{equation*}
|
||
|
|
K(x,y) = K_1(x,y) = c \frac{(x-y)^\perp}{|x-y|^2} .
|
||
|
|
\end{equation*}
|
||
|
|
Here $x^\perp = (-x_2,x_1)$, and $c$ is a constant, probably something
|
||
|
|
like $1/2\pi$. In any case we will set it to be one, which in effect
|
||
|
|
is rescaling the time variable, so we don't need to worry about it.
|
||
|
|
|
||
|
|
We can use this as a basis to find $K$ on the unit disk
|
||
|
|
$\Omega = \Delta = \{x:|x|<1\}$. It turns out to be
|
||
|
|
\begin{equation*}
|
||
|
|
K_2(x,y) = K_1(x,y) - K_1(x,y^*) ,
|
||
|
|
\end{equation*}
|
||
|
|
where $y^* = y/|y|^2$ is called the reflection of $y$ about the
|
||
|
|
boundary of the unit disk.
|
||
|
|
|
||
|
|
Another example is if we have a bijective analytic function
|
||
|
|
$p:\Delta \to {\mathbb C}$, and we let $\Omega = p(\Delta)$.
|
||
|
|
(Here we think of $\Delta$ as a subset of $\mathbb C$, that is,
|
||
|
|
we are identifying the plane with the set of complex numbers.)
|
||
|
|
In that case we get
|
||
|
|
\[ K_p(p(x),p(y)) = K_2(x,y)/|p'(x)|^2 .\]
|
||
|
|
Our simulation considers the last case. Examples of such
|
||
|
|
analytic functions include series
|
||
|
|
$p(x) = x + \sum_{n=2}^\infty c_n x^n$, where
|
||
|
|
$\sum_{n=2}^\infty n |c_n| \le 1$.
|
||
|
|
(Thanks to David Ullrich for pointing this out to me.)
|
||
|
|
|
||
|
|
\section{The Simulation}
|
||
|
|
|
||
|
|
Now let's get to decribing the simulation. We assume a rather
|
||
|
|
unusual initial distribution for the vorticity --- that the
|
||
|
|
vorticity is a finite sum of dirac delta masses.
|
||
|
|
\[ w(0,x) = \sum_{k=1}^N w_k \delta(x-x_k(0)) .\]
|
||
|
|
Here $x_k(0)$ is the initial place where the points
|
||
|
|
of vorticity are concentrated, with values $w_k$.
|
||
|
|
Then at time $t$, the vorticity becomes
|
||
|
|
\[ w(t,x) = \sum_{k=1}^N w_k \delta(x-x_k(t)) .\]
|
||
|
|
The points of fluid $x_k(t)$ are pushed by the
|
||
|
|
flow, that is, $x_k(t) = \varphi(t,x_k(0))$, or
|
||
|
|
\[ \frac{\partial}{\partial t} x_k(t) = u(t,x_k(t)) .\]
|
||
|
|
Putting this all together, we finally obtain the equations
|
||
|
|
\[ \frac{\partial}{\partial t} x_k = \alpha_k \]
|
||
|
|
where
|
||
|
|
\[ \alpha_k = \sum_{l=1}^N w_l K(x_k,x_l) .\]
|
||
|
|
This is the equation that our simulation solves.
|
||
|
|
|
||
|
|
In fact, in our case, where the domain is $p(\Delta)$,
|
||
|
|
the points are described by points
|
||
|
|
$\tilde x_k$, where $x_k = p(\tilde x_k)$. Then
|
||
|
|
the equations become
|
||
|
|
\begin{eqnarray}
|
||
|
|
\label{tildex-p1}
|
||
|
|
\frac{\partial}{\partial t} \tilde x_k &=& \tilde\alpha_k \\
|
||
|
|
\label{tildex-p2}
|
||
|
|
\tilde\alpha_k &=& \frac1{|p'(\tilde x_k)|^2}
|
||
|
|
\sum_{l=1}^N w_l K_2(\tilde x_k,\tilde x_l) .
|
||
|
|
\end{eqnarray}
|
||
|
|
|
||
|
|
We solve this $2N$ system of equations using standard
|
||
|
|
numerical methods, in our case, using the second order midpoint method
|
||
|
|
for the first step, and thereafter using the second order Adams-Bashforth
|
||
|
|
method. (See for example the book
|
||
|
|
by Burden and Faires \cite{BF}).
|
||
|
|
|
||
|
|
\section{The Program - Data Structures}
|
||
|
|
|
||
|
|
The computer program solves equation (\ref{tildex-p1}), and displays
|
||
|
|
the results on the screen, with a boundary. All the information
|
||
|
|
for solving the equation and displaying the output is countained
|
||
|
|
in the structure {\tt euler2dstruct}. Let us describe some of
|
||
|
|
the fields in {\tt euler2dstruct}.
|
||
|
|
The points $\tilde x_k$ are contained
|
||
|
|
in {\tt double *x}: with the coordinates of
|
||
|
|
$\tilde x_k$ being the two numbers
|
||
|
|
{\tt x[2*k+0]}, {\tt x[2*k+1]}. The values $w_k$ are contained
|
||
|
|
in {\tt double *w}. The total number of points is
|
||
|
|
{\tt int N}. (But only the first {\tt int Nvortex} points
|
||
|
|
have $w_k \ne 0$.) The coefficients of the analytic function
|
||
|
|
(in our case a polynomial) $p$
|
||
|
|
are contained in {\tt double p\_coef[2*(deg\_p-1)]} --- here
|
||
|
|
{\tt deg\_p} is the degree of $p$, and the real and imaginary
|
||
|
|
parts of the coefficient
|
||
|
|
$c_n$ is contained in {\tt p\_coef[2*(n-2)+0]} and {\tt p\_coef[2*(n-2)+1]}.
|
||
|
|
|
||
|
|
\section{Data Initialization}
|
||
|
|
|
||
|
|
The program starts in the function {\tt init\_euler2d}. After allocating
|
||
|
|
the memory for the data, and initialising some of the temporary variables
|
||
|
|
required for the numerical solving program, it randomly assigns the
|
||
|
|
coefficients of $p$, making sure that $\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$.
|
||
|
|
Then the program figures out how to draw the boundary, and what rescaling
|
||
|
|
of the data is required to draw it on the screen. (This uses the
|
||
|
|
function {\tt calc\_p} which calculates $p(x)$.)
|
||
|
|
|
||
|
|
Next, it randomly assigns the initial values of $\tilde x_k$. We want
|
||
|
|
to do this in such a way so that the points are uniformly spread over the
|
||
|
|
domain. Let us first consider the case when the domain is the unit circle
|
||
|
|
$\Delta$. In that case the proportion of points that we would expect
|
||
|
|
inside the circle of radius $r$ would be proportional to $r^2$. So
|
||
|
|
we do it as follows:
|
||
|
|
\[ r = \sqrt{R_{0,1}},\quad \theta = R_{-\pi,\pi}, \quad
|
||
|
|
\tilde x_k = r (\cos \theta, \sin \theta) .\]
|
||
|
|
Here, and in the rest of this discussion, $R_{a,b}$ is a function
|
||
|
|
that returns a random variable uniformly distributed over the interval
|
||
|
|
$[a,b]$.
|
||
|
|
|
||
|
|
This works fine for $\Delta$, but for $p(\Delta)$, the points
|
||
|
|
$p(\tilde x_k)$ are not uniformly distributed over $p(\Delta)$,
|
||
|
|
but are distributed with a density proportional to
|
||
|
|
$1/|p'(\tilde x_k)|^2$. So to restore the uniform density we need
|
||
|
|
to reject this value of $\tilde x_k$ with probability proportional
|
||
|
|
to $|p'(\tilde x_k)|^2$. Noticing that the condition
|
||
|
|
$\sum_{n=2}^{\tt deg\_p} n |c_n| = 1$ implies that
|
||
|
|
$|p'(\tilde x_k)| \le 2$, we
|
||
|
|
do this by rejecting if $|p'(\tilde x_k)|^2 < R_{0,4}$.
|
||
|
|
(This makes use of the function {\tt calc\_mod\_dp2} which calculates
|
||
|
|
$|p'(x)|^2$.)
|
||
|
|
|
||
|
|
\section{Solving the Equation}
|
||
|
|
|
||
|
|
The main loop of the program is in the function {\tt draw\_euler2d}.
|
||
|
|
Most of the drawing operations are contained in this function, and
|
||
|
|
the numerical aspects are sent to the function {\tt ode\_solve}.
|
||
|
|
But there is an aspect of this that I would like
|
||
|
|
to discuss in the next section, and so we will look at a simple method for
|
||
|
|
numerically solving differential equations.
|
||
|
|
|
||
|
|
The Euler Method
|
||
|
|
(nothing to do with the Euler Equation), is as
|
||
|
|
follows. Pick a small number $h$ --- the time step (in
|
||
|
|
the program call {\tt delta\_t}). Then we approximate
|
||
|
|
the solution of the equation:
|
||
|
|
\begin{equation}
|
||
|
|
\label{method-simple}
|
||
|
|
\tilde x_k(t+h) = \tilde x_k(t) + h \tilde\alpha_k(t) .
|
||
|
|
\end{equation}
|
||
|
|
The more sophisticated methods we use are variations of
|
||
|
|
the Euler Method, and so the discussion in the following section
|
||
|
|
still applies.
|
||
|
|
|
||
|
|
In the program, the quantities $\tilde\alpha_k$, given by
|
||
|
|
equations (\ref{tildex-p2}) are calculated by the function
|
||
|
|
{\tt derivs}
|
||
|
|
(which in turns calls {\tt calc\_all\_mod\_dp2} to
|
||
|
|
calculate $|p'(\tilde x_k)|^2$ at all the points).
|
||
|
|
|
||
|
|
|
||
|
|
\section{Subtle Perturbation}
|
||
|
|
|
||
|
|
Added later: the scheme described here seems to not be that effective,
|
||
|
|
so now it is not used.
|
||
|
|
|
||
|
|
One problem using a numerical scheme such as the Euler Method occurs
|
||
|
|
when the points $\tilde x_k$ get close to the boundary
|
||
|
|
of $\Delta$. In that case, it is possible that the new
|
||
|
|
points will be pushed outside of the boundary. Even if they
|
||
|
|
are not pushed out of the boundary, they may be much closer
|
||
|
|
or farther from the boundary than they should be.
|
||
|
|
Our system of equations is very sensitive to how close points
|
||
|
|
are to the boundary --- points with non-zero vorticity
|
||
|
|
(``vortex points'') that are close to the boundary travel
|
||
|
|
at great speed alongside the boundary, with speed that is
|
||
|
|
inversely proportional to the distance from the boundary.
|
||
|
|
|
||
|
|
A way to try to mitigate this problem is something that I call
|
||
|
|
``subtle perturbation.''
|
||
|
|
We map the points in
|
||
|
|
the unit disk to points in the plane using the map
|
||
|
|
\begin{equation*}
|
||
|
|
F(x) = f(|x|) \frac x{|x|} ,
|
||
|
|
\end{equation*}
|
||
|
|
where $f:[0,1]\to[0,\infty]$ is an increasing continuous
|
||
|
|
bijection. It turns out that a good choice is
|
||
|
|
\begin{equation*}
|
||
|
|
f(t) = -\log(1-t) .
|
||
|
|
\end{equation*}
|
||
|
|
(The reason for this is that points close to each other
|
||
|
|
that are a distance
|
||
|
|
about $r$ from the boundary will be pushed around so that
|
||
|
|
their distance from each other is about multiplied by the
|
||
|
|
derivative of $\log r$, that is, $1/r$.)
|
||
|
|
Note that the inverse of this function is given by
|
||
|
|
\begin{equation*}
|
||
|
|
F^{-1}(x) = f^{-1}(|x|) \frac x{|x|} ,
|
||
|
|
\end{equation*}
|
||
|
|
where
|
||
|
|
\begin{equation*}
|
||
|
|
f^{-1}(t) = 1-e^{-t} .
|
||
|
|
\end{equation*}
|
||
|
|
|
||
|
|
So what we could do is the following: instead of working with
|
||
|
|
the points $\tilde x_k$, we could work instead with the points
|
||
|
|
$y_k = F(\tilde x_k)$. In effect this is what we do.
|
||
|
|
Instead of performing the computation (\ref{method-simple}),
|
||
|
|
we do the calculation
|
||
|
|
\begin{equation*}
|
||
|
|
y_k = F(\tilde x_k(t)) + h {\cal A}(\tilde x_k) \tilde\alpha_k(t)
|
||
|
|
\end{equation*}
|
||
|
|
where
|
||
|
|
${\cal A}(x)$ is the matrix of partial derivatives of $F$:
|
||
|
|
\begin{equation*}
|
||
|
|
{\cal A}(x) =
|
||
|
|
\frac{f(|x|)}{|x|}
|
||
|
|
\left[
|
||
|
|
\begin{matrix}
|
||
|
|
1 & 0\\
|
||
|
|
0 & 1
|
||
|
|
\end{matrix}
|
||
|
|
\right]
|
||
|
|
+ \frac1{|x|}
|
||
|
|
\left(\frac{f'(|x|)}{|x|} - \frac{f(|x|)}{|x|^2}\right)
|
||
|
|
\left[
|
||
|
|
\begin{matrix}
|
||
|
|
x_{1}^2 & x_{1} x_{2}\\
|
||
|
|
x_{1} x_{2} & x_{2}^2
|
||
|
|
\end{matrix}
|
||
|
|
\right],
|
||
|
|
\end{equation*}
|
||
|
|
and then compute
|
||
|
|
\begin{equation*}
|
||
|
|
\tilde x_k(t+h) = F^{-1}(y_k).
|
||
|
|
\end{equation*}
|
||
|
|
These calculations are done in the function {\tt perturb}, if
|
||
|
|
the quantity {\tt SUBTLE\_PERTURB} is set.
|
||
|
|
|
||
|
|
\section{Drawing the Points}
|
||
|
|
|
||
|
|
As we stated earlier, most of the drawing functions are contained
|
||
|
|
in the function {\tt draw\_euler2d}. If the variable
|
||
|
|
{\tt hide\_vortex} is set (and the function {\tt init\_euler2d}
|
||
|
|
will set this with probability $3/4$), then we only display
|
||
|
|
the points $\tilde x_k$ for ${\tt Nvortex} < k \le N$. If
|
||
|
|
{\tt hide\_vortex} is not set, then the ``vortex points''
|
||
|
|
$\tilde x_k$ ($1 \le k \le {\tt Nvortex}$) are displayed in white.
|
||
|
|
In fact the points $p(\tilde x_k)$ are what are put onto the screen,
|
||
|
|
and for this we make use of the function {\tt calc\_all\_p}.
|
||
|
|
|
||
|
|
\section{Addition to Program: Changing the Power Law}
|
||
|
|
|
||
|
|
A later addition to the program adds an option {\tt eulerpower},
|
||
|
|
which allows one to change the power law that describes how
|
||
|
|
the vortex points influence other points. In effect, if this
|
||
|
|
option is set with the value $m$, then the Biot-Savart Kernel
|
||
|
|
is replace by
|
||
|
|
$$ K_1(x,y) = \frac{(x-y)^\perp}{|x-y|^{m+1}}, $$
|
||
|
|
and
|
||
|
|
$$ K_2(x,y) = K_1(x,y) - |y|^{1-m} K_1(x,y) .$$
|
||
|
|
So for example, setting $m=2$ corresponds to the
|
||
|
|
quasi-geostrophic equation. (I haven't yet figured out
|
||
|
|
what $K_p$ should be, so if $m \ne 1$ we use the unit circle
|
||
|
|
as the boundary.)
|
||
|
|
|
||
|
|
\begin{thebibliography}{9}
|
||
|
|
|
||
|
|
\bibitem{BF} Richard L. Burden, J. Douglas Faires, Numerical Analysis,
|
||
|
|
sixth edition, Brooks/Cole, 1996.
|
||
|
|
|
||
|
|
\bibitem{C} Alexandre J. Chorin, Vorticity and Turbulence,
|
||
|
|
Applied Mathematical Sciences, Vol 103, Springer Verlag, 1994.
|
||
|
|
|
||
|
|
\end{thebibliography}
|
||
|
|
|
||
|
|
\end{document}
|