<br>
-<i>This program was contributed by Luca Heltai.
+<i>This program was contributed by Luca Heltai.
</i>
<a name="Intro"></a>
The three dimensional solution generated by this example program is
wrong. There is a factor of $1/2$ that appears, and we haven't figured
out where it comes from yet. For the moment this issue is "solved" by
-multiplying the wind function in the parameter file by two.
+multiplying the wind function in the parameter file by two.
If you think you have spotted the mistake, please let us know.
\f{align*}
\frac{\partial }{\partial t}\mathbf{v} + (\mathbf{v}\cdot\nabla)\mathbf{v}
&=
- -\frac{1}{\rho}\nabla p + \mathbf{g}
+ -\frac{1}{\rho}\nabla p + \mathbf{g}
\qquad &\text{in } \mathbb{R}^n \backslash \Omega
\\
- \nabla \cdot \mathbf{v}&=0
+ \nabla \cdot \mathbf{v}&=0
&\text{in } \mathbb{R}^n\backslash\Omega
\f}
where the fluid density $\rho$ and the acceleration $\mathbf{g}$ due
to external forces are given and the velocity $\mathbf{v}$ and the
pressure $p$ are the unknowns. Here $\Omega$ is a closed bounded
-region representing the body around which the fluid moves.
+region representing the body around which the fluid moves.
The above equations can be derived from Navier-Stokes equations
assuming that the effects due to viscosity are negligible compared to
-\frac{1}{\rho}\nabla p
\qquad &\text{in } \mathbb{R}^n \backslash \Omega
\\
- \nabla \cdot \mathbf{v}&=0
+ \nabla \cdot \mathbf{v}&=0
&\text{in } \mathbb{R}^n\backslash\Omega
\f}
coordinate system moves along with the body whereas the background
fluid is at rest at infinity. Notice that we define the normal
$\mathbf{n}$ as the <i>outer</i> normal to the domain $\Omega$, which
-is the opposite of the outer normal to the integration domain.
+is the opposite of the outer normal to the integration domain.
For both stationary and non stationary flow, the solution process
starts by solving for the velocity in the second equation and
substituting in the first equation in order to find the pressure.
The solution of the stationary Euler equations is typically performed
in order to understand the behavior of the given (possibly complex)
-geometry when a prescribed motion is enforced on the system.
+geometry when a prescribed motion is enforced on the system.
The first step in this process is to change the frame of reference from a
coordinate system moving along with the body to one in which the body moves
\label{laplace}
\Delta\phi &= 0 \qquad &&\text{in}\ \mathbb{R}^n\backslash\Omega,
\\
- \mathbf{n}\cdot\nabla\phi &= -\mathbf{n}\cdot\mathbf{v}_\infty
+ \mathbf{n}\cdot\nabla\phi &= -\mathbf{n}\cdot\mathbf{v}_\infty
&& \text{on}\ \partial\Omega
\f}
while the momentum equation reduces to Bernoulli's equation that expresses the
So we can solve the problem by solving the Laplace equation for the
potential. We recall that the following functions, called fundamental
-solutions of the Laplace equation,
+solutions of the Laplace equation,
-\f[ \begin{aligned}
+\f[ \begin{aligned}
\label{eq:3} G(\mathbf{y}-\mathbf{x}) = &
-\frac{1}{2\pi}\ln|\mathbf{y}-\mathbf{x}| \qquad && \text{for } n=2 \\
G(\mathbf{y}-\mathbf{x}) = &
\frac{1}{4\pi}\frac{1}{|\mathbf{y}-\mathbf{x}|}&& \text{for } n=3,
-\end{aligned}
-\f]
+\end{aligned}
+\f]
-satisfy in a distributional sense the equation:
+satisfy in a distributional sense the equation:
-\f[
--\Delta_y G(\mathbf{y}-\mathbf{x}) = \delta(\mathbf{y}-\mathbf{x}),
-\f]
+\f[
+-\Delta_y G(\mathbf{y}-\mathbf{x}) = \delta(\mathbf{y}-\mathbf{x}),
+\f]
where the derivative is done in the variable $\mathbf{y}$. By using
the usual Green identities, our problem can be written on the boundary
\f[\label{green}
\int_{\omega}
(-\Delta u)v\,dx + \int_{\partial\omega} \frac{\partial u}{\partial \mathbf{n} }v \,ds
- =
+ =
\int_{\omega}
- (-\Delta v)u\,dx + \int_{\partial\omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
+ (-\Delta v)u\,dx + \int_{\partial\omega} u\frac{\partial v}{\partial \mathbf{n}} \,ds,
\f]
where $\mathbf{n}$ is the normal to the surface of $\omega$ pointing
respectively, as long as $\mathbf{x}$ is chosen in the region
$\mathbb{R}^n\backslash\Omega$, we obtain:
\f[
- \phi(\mathbf{x}) -
- \int_{\Gamma\cup\Gamma_\infty}\frac{\partial G(\mathbf{y}-\mathbf{x})}{\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y
+ \phi(\mathbf{x}) -
+ \int_{\Gamma\cup\Gamma_\infty}\frac{\partial G(\mathbf{y}-\mathbf{x})}{\partial \mathbf{n}_y}\phi(\mathbf{y})\,ds_y
=
-\int_{\Gamma\cup\Gamma_\infty}G(\mathbf{y}-\mathbf{x})\frac{\partial \phi}{\partial \mathbf{n}_y}(\mathbf{y})\,ds_y
\qquad \forall\mathbf{x}\in \mathbb{R}^n\backslash\Omega
\f[
-\int_{\Gamma_\infty} \frac{\partial G(\mathbf{y}-\mathbf{x})}
-{\partial \mathbf{n}_y}\phi_\infty \,ds_y =
-\lim_{r\to\infty} \int_{\partial B_r(0)} \frac{\mathbf{r}}{r} \cdot \nabla G(\mathbf{y}-\mathbf{x})
+{\partial \mathbf{n}_y}\phi_\infty \,ds_y =
+\lim_{r\to\infty} \int_{\partial B_r(0)} \frac{\mathbf{r}}{r} \cdot \nabla G(\mathbf{y}-\mathbf{x})
\phi_\infty \,ds_y = -\phi_\infty.
\f]
$\mathbf{n}\cdot\nabla\phi = -\mathbf{n}\cdot\mathbf{v}_\infty$.
Consequently,
\f[
- \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
+ \phi(\mathbf{x}) - (D\phi)(\mathbf{x}) = \phi_\infty +
\left(S[\mathbf{n}\cdot\mathbf{v}_\infty]\right)(\mathbf{x})
\qquad \forall\mathbf{x} \in \mathbb{R}^n\backslash\Omega.
\f]
at points where the boundary has a corner or an edge.
Substituting the single and double layer operators we get:
-\f[
- \alpha(\mathbf{x}) \phi(\mathbf{x})
+\f[
+ \alpha(\mathbf{x}) \phi(\mathbf{x})
+ \frac{1}{2\pi}\int_{\partial \Omega} \frac{
(\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }{ |\mathbf{y}-\mathbf{x}|^2 }
- \phi(\mathbf{x}) \,ds_y
+ \phi(\mathbf{y}) \,ds_y
= \phi_\infty
-\frac{1}{2\pi}\int_{\partial \Omega} \ln|\mathbf{y}-\mathbf{x}| \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
-\f]
+\f]
for two dimensional flows and
-\f[
- \alpha(\mathbf{x}) \phi(\mathbf{x})
+\f[
+ \alpha(\mathbf{x}) \phi(\mathbf{x})
+ \frac{1}{4\pi}\int_{\partial \Omega} \frac{ (\mathbf{y}-\mathbf{x})\cdot\mathbf{n}_y }{ |\mathbf{y}-\mathbf{x}|^3 }\phi(\mathbf{y})\,ds_y
= \phi_\infty +
\frac{1}{4\pi}\int_{\partial \Omega} \frac{1}{|\mathbf{y}-\mathbf{x}|} \, \mathbf{n}\cdot\mathbf{v_\infty}\,ds_y
-\f]
+\f]
for three dimensional flows, where the normal derivatives of the fundamental
solutions have been written in a form that makes computation easier. In either
case, $\phi$ is the solution of an integral equation posed entirely on the
can compute (analytically, if we want) the solution $\phi(\mathbf{x})$ in all
of $\mathbb{R}^n\backslash\Omega$. To this end, recall that we had
\f[
- \phi(\mathbf{x})
+ \phi(\mathbf{x})
=
\phi_\infty +
(D\phi)(\mathbf{x})
where now we have everything that is on the right hand side ($S$ and $D$ are
integrals we can evaluate, the normal velocity on the boundary is given, and
$\phi$ on the boundary we have just computed). Finally, we can then recover
-the velocity as $\mathbf{\tilde v}=\nabla \phi$.
+the velocity as $\mathbf{\tilde v}=\nabla \phi$.
Notice that the evaluation of the above formula for $\mathbf{x} \in
\Omega$ should yield zero as a result, since the integration of the
&=
\partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
-
- \partial_j [(\partial_j\phi+v_{\infty,j})] \partial_i\phi
+ \partial_j [(\partial_j\phi+v_{\infty,j})] \partial_i\phi
\\
&=
\partial_j [(\partial_j\phi+v_{\infty,j}) \partial_i\phi]
@f}
We now only need to massage that last term a bit more. Using the product rule,
-we get
+we get
@f{align*}
\partial_j [\partial_i (\partial_j\phi) \phi]
&=
\partial_j [\partial_i (\partial_j\phi) \phi]
\\
&=
- \partial_i
+ \partial_i
\left[
(\partial_j\phi)(\partial_j \phi) + v_{\infty,j} \partial_j \phi
- -
+ -
\frac 12 (\partial_j\phi)(\partial_j\phi)
\right],
\\
&=
- \partial_i
+ \partial_i
\left[
\frac 12 (\partial_j\phi)(\partial_j \phi) + v_{\infty,j} \partial_j \phi
\right],
@f[
\mathbf{v}\cdot\nabla\mathbf{v}
=
- \nabla
+ \nabla
\left[
\frac 12 \mathbf{\tilde v}^2
+ \mathbf{v}_{\infty} \cdot \mathbf{\tilde v}
@f[
p
=
- -\rho
+ -\rho
\left[
\frac 12 \mathbf{\tilde v}^2
+ \mathbf{v}_{\infty} \cdot \mathbf{\tilde v}
\right]
=
- -\rho
+ -\rho
\left[
\frac 12 \mathbf{v}^2
-
quadrilaterals if $n=3$. We will call each individual segment or
quadrilateral an <i>element</i> or <i>cell</i>, independently of the
dimension $n$ of the surrounding space $\mathbb{R}^n$.
-We define the finite dimensional space $V_h$ as
+We define the finite dimensional space $V_h$ as
\f[
\label{eq:definition-Vh}
- V_h := \{ v \in C^0(\Gamma) \text{ s.t. } v|_{K_i} \in \mathcal{Q}^1(K_i),
+ V_h := \{ v \in C^0(\Gamma) \text{ s.t. } v|_{K_i} \in \mathcal{Q}^1(K_i),
\forall i\},
\f]
with basis functions $\psi_i(\mathbf{x})$ for which we will use the usual FE_Q
finite element, with the catch that this time it is defined on a manifold of
codimension one (which we do by using the second template argument that is
usually defaulted to equal the first; here, we will create objects
-<code>FE_Q@<dim+1,dim@></code> to indicate that we have <code>dim-1</code>
+<code>FE_Q@<dim-1,dim@></code> to indicate that we have <code>dim-1</code>
dimensional cells in a <code>dim</code> dimensional space).
-An element $\phi_h$ of $V_h$ is uniquely
+An element $\phi_h$ of $V_h$ is uniquely
identified by the vector $\boldsymbol{\phi}$ of its coefficients
$\phi_i$, that is:
\f[
\label{eq:definition-of-element}
- \phi_h(\mathbf{x}) := \phi_i \psi_i(\mathbf{x}), \qquad
+ \phi_h(\mathbf{x}) := \phi_i \psi_i(\mathbf{x}), \qquad
\boldsymbol{\phi} := \{ \phi_i \},
\f]
where summation is implied over repeated indexes. Note that we could use
discontinuous elements here — in fact, there is no real reason to use
continuous ones since the integral formulation does not
-imply any derivatives on our trial functions so continuity is unnecessary,
+imply any derivatives on our trial functions so continuity is unnecessary,
and often in the literature only piecewise constant elements are used.
<h3> Collocation boundary element method </h3>
such that the following $n\_dofs$ equations are satisfied:
\f{align*}
- \alpha(\mathbf{x}_i) \phi_h(\mathbf{x}_i)
+ \alpha(\mathbf{x}_i) \phi_h(\mathbf{x}_i)
- \int_{\Gamma_y} \frac{ \partial G(\mathbf{y}-\mathbf{x}_i)}{\partial\mathbf{n}_y }
- \phi_h(\mathbf{y}) \,ds_y =
- \int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \,
- \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y
+ \phi_h(\mathbf{y}) \,ds_y =
+ \int_{\Gamma_y} G(\mathbf{y}-\mathbf{x}_i) \,
+ \mathbf{n}_y\cdot\mathbf{v_\infty} \,ds_y
,
\f}
\f[
\begin{aligned}
-\mathbf{A}_{ij}&=
+\mathbf{A}_{ij}&=
\alpha(\mathbf{x}_i) \psi_j(\mathbf{x}_i)
-= 1+\int_\Gamma
-\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
+= 1+\int_\Gamma
+\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
\psi_j(\mathbf{x}_i)
\\
\mathbf{N}_{ij}&= - \int_\Gamma
\frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}
- \psi_j(\mathbf{y}) \,ds_y
+ \psi_j(\mathbf{y}) \,ds_y
\\
-\mathbf{b}_i&= \int_\Gamma
+\mathbf{b}_i&= \int_\Gamma
G(\mathbf{y}-\mathbf{x}_i) \, \mathbf{n}_y\cdot\mathbf{v_\infty}
ds_y.
\end{aligned}
$\mathbf{A}$ is diagonal with entries
\f[
\mathbf{A}_{ii}
- =
- 1+\int_\Gamma
- \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
+ =
+ 1+\int_\Gamma
+ \frac{\partial G(\mathbf{y}-\mathbf{x}_i)}{\partial \mathbf{n}_y}\,ds_y
=
1-\sum_j N_{ij},
\f]
where we have used that $\sum_j \psi_j(\mathbf{y})=1$ for the usual Lagrange
-elements.
+elements.
With this choice of collocation points, the computation of the entries
of the matrices $\mathbf{A}$, $\mathbf{N}$ and of the right hand side
$\mathbf{b}$ requires the evaluation of singular integrals on the
bi-linear (in three dimensions) transformation of the reference
boundary element $\hat K := [0,1]^{n-1}$, and we perform the integrations after a
change of variables from the real element $K_i$ to the reference
-element $\hat K$.
+element $\hat K$.
<h3> Treating the singular integrals. </h3>
goes to zero when $\mathbf{x}=\mathbf{y}$, the numerator is always
zero because $\mathbf{n}_y$ and $(\mathbf{y}-\mathbf{x})$ are
orthogonal (on our polygonal approximation of the boundary of $\Omega$), and
-the only singular integral arises in the computation
+the only singular integral arises in the computation
of $\mathbf{b}_i$ on the i-th element of $\mathcal{T}_h$:
\f[
\frac{1}{\pi}
significantly sparser representation of these matrices that also facilitates
rapid evaluations of the scalar product between vectors and matrices. This not
being the goal of this program, we leave this for more sophisticated
-implementations.
+implementations.
<h3>Implementation</h3>
However, this doesn't have to be so: in the current example, we will for
example want to solve on the surface of a sphere, which is a two-dimensional
manifold embedded in a three-dimensional space. Consequently, the right class
-will be <code>Triangulation@<2,3@></code>, and correspondingly we will use
+will be <code>Triangulation@<2,3@></code>, and correspondingly we will use
<code>DoFHandler@<2,3@></code> as the DoF handler class and
<code>FE_Q@<2,3@></code> for finite elements.
projects / dealii.git / commitdiff