Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
SWIFTsim
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Deploy
Releases
Model registry
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
SWIFT
SWIFTsim
Merge requests
!494
Add EOM derivation
Code
Review changes
Check out branch
Download
Patches
Plain diff
Merged
Add EOM derivation
JB_update_sph_documentation
into
master
Overview
7
Commits
3
Changes
3
Merged
Josh Borrow
requested to merge
JB_update_sph_documentation
into
master
7 years ago
Overview
7
Commits
3
Changes
3
Expand
In this merge request I add:
Derivation of the Lagrange multipliers
Derivation of the SPH EOM
Derivation of the 'f_{ij}' factors.
0
0
Merge request reports
Compare
master
version 1
eda861ae
7 years ago
master (base)
and
latest version
latest version
82c644a0
3 commits,
7 years ago
version 1
eda861ae
2 commits,
7 years ago
3 files
+
230
−
8
Inline
Compare changes
Side-by-side
Inline
Show whitespace changes
Show one file at a time
Files
3
Search (e.g. *.vue) (Ctrl+P)
theory/SPH/Derivation/sph_derivation.tex
0 → 100644
+
219
−
0
Options
Following
\citet
{
hopkins2013
}
we use the Lagrangian formalism to determine the
equation of motion for the SPH particles. The derivation found in
\citet
{
hopkins2013
}
is somewhat sparse, and so we reproduce it here with more
steps for clarity.
The following derivation is underpinned by the idea of there being two
independent ways of defining the volume associated with a particle in SPH. The
first is the volume associated with the thermodynamical system (
$
\Delta
V
$
),
from the first law of thermodynamics, and the second being the volume around the
particle in which we conserve an effective neighbour number (
$
\Delta
\tilde
{
V
}$
).
These two need not necessarily be linked in any way.
We begin with the SPH lagragian,
\begin{align}
L(q,
\dot
{
q
}
) =
\frac
{
1
}{
2
}
\sum
^
N
_{
i=1
}
m
_
i
\dot
{
r
}^
2
_
i
-
\sum
^
N
_{
i=1
}
m
_
i u
_
i,
\label
{
eqn:sph:derivation:sphlagrangian
}
\end{align}
and the first law of thermodynamics,
\begin{align}
\left
.
\frac
{
\partial
u
_
i
}{
\partial
q
_
i
}
\right
|
_
A =
-
\frac
{
P
_
i
}{
m
_
i
}
\frac
{
\partial
\Delta
V
_
i
}{
\partial
q
_
i
}
,
\label
{
eqn:sph:derivation:firstlaw
}
\end{align}
where
$
\mathbf
{
q
}
=
(
\mathbf
{
r
}_
1
, ...,
\mathbf
{
r
}_
N, h
_
i, ..., h
_
N
)
$
are the
generalised coordinates of the particles, and the derivative of the internal
energy
$
u
_
i
$
is taken at fixed entropy
$
A
$
. This gives the first concept of
`volume'
$
\Delta
V
_
i
$
that the particles occupy.
As mentioned earlier, the particles also have a volume associated with the
spread of their neighbours and hence their smoothing length. We can write this as
a constraint equation,
\begin{align}
\phi
_
i(
\mathbf
{
q
}
) =
\kappa
h
_
i
^{
n
_
d
}
\frac
{
1
}{
\Delta
\tilde
{
V
}}
- N
_{
ngb
}
= 0,
\label
{
eqn:sph:derivation:constraint
}
\end{align}
where
$
N
_{
ngb
}$
is the effecitve neighbour number,
$
\kappa
$
is the volume of
the unit sphere (
$
\kappa
_{
3
D
}
=
4
\pi
/
3
$
), and
$
n
_
d
$
is the number of spatial
dimensions considered in the problem. It is important to note that
$
N
_{
ngb
}$
need not be an integer quantity.
\subsection
{
Lagrange Multipliers
}
With these components in hand, we can use the technique of Lagrange multipliers
to enforce the constraint equation, with
\begin{align}
\frac
{
\mathrm
{
d
}}{
\mathrm
{
d
}
t
}
\frac
{
\partial
L
}{
\partial
\dot
{
q
}_
i
}
-
\frac
{
\partial
L
}{
\partial
q
_
i
}
=
\sum
^{
N
}_{
j=1
}
\lambda
_
j
\frac
{
\partial
\phi
_
i
}{
\partial
q
_
j
}
,
\label
{
eqn:sph:derivation:lmsum
}
\end{align}
where the
$
\lambda
_
j
$
are the lagrange multipliers. We use the second half of
these equations (i.e.
$
q
_
i
=
h
_
i
$
) to constrain
$
\lambda
_
j
$
. The differentials
with respect to the smoothing lengths:
\begin{align}
\frac
{
\partial
L
}{
\partial
\dot
{
h
}_
i
}
= 0,
\quad
\frac
{
\partial
L
}{
\partial
h
_
i
}
=
-
\sum
^
N
_{
j=1
}
m
_
j
\frac
{
\partial
u
_
j
}{
\partial
h
_
i
}
=
-m
_
i
\frac
{
\partial
u
_
i
}{
\partial
h
_
i
}
.
\end{align}
As all terms where
$
i
\neq
j
=
0
$
. For the constraint equation we find
\begin{align}
\frac
{
\partial
\phi
_
j
}{
\partial
h
_
i
}
=
\kappa
n
_
d h
_
j
^{
n
_
d -1
}
\frac
{
\partial
h
_
j
}{
\partial
h
_
i
}
\frac
{
1
}{
\Delta
\tilde
{
V
}_
j
}
+
\kappa
h
_
j
^{
n
_
d
}
\frac
{
\partial
\Delta
\tilde
{
V
}_
j
}{
\partial
h
_
i
}
\frac
{
1
}{
\Delta
\tilde
{
V
}_
j
^
2
}
,
\end{align}
which clearly reduces to
\begin{align}
\frac
{
\partial
\phi
_
j
}{
\partial
h
_
i
}
=
\kappa
\left
(n
_
d h
_
j
^{
n
_
d - 1
}
\frac
{
1
}{
\Delta
\tilde
{
V
_
j
}}
+
h
_
j
^{
n
_
d
}
\frac
{
\partial
\Delta
\tilde
{
V
_
j
}}{
\partial
h
_
i
}
\frac
{
1
}{
\Delta
\tilde
{
V
_
j
}^
2
}
\right
)
\delta
_{
ij
}
,
\end{align}
The first law of thermodynamics gives us an expression for
${
\partial
u
_
i
}
/
{
\partial
h
_
i
}$
,
\begin{align}
\frac
{
\partial
u
_
i
}{
\partial
h
_
i
}
=
-
\frac
{
P
_
i
}{
m
_
i
}
\frac
{
\partial
\Delta
V
_
i
}{
\partial
h
_
i
}
.
\end{align}
Putting all this together we find that
\begin{align}
P
_
i
\frac
{
\partial
\Delta
V
_
i
}{
\partial
h
_
i
}
=
\sum
^
N
_{
j=1
}
\kappa
\lambda
_
j
&
\left
( n
_
d h
_
j
^{
n
_
d -1
}
\frac
{
1
}{
\Delta
\tilde
{
V
_
j
}}
\right
.
\nonumber\\
&
\left
. + ~ h
_
j
^{
n
_
d
}
\frac
{
\partial
\Delta
\tilde
{
V
_
j
}}{
\partial
h
_
i
}
\frac
{
1
}{
\Delta
\tilde
{
V
_
j
}^
2
}
\right
)
\delta
_{
ij
}
.
\end{align}
This expression can now the rearranged for
$
\lambda
_
i
$
,
\begin{align}
\lambda
_
i =
\left
(
\frac
{
P
_
i
}{
\kappa
}
\frac
{
\Delta
\tilde
{
V
}^
2
_
i
}{
h
^{
n
_
d
}_
i
}
\right
)
\frac
{
h
_
i
}{
n
_
d
\Delta
\tilde
{
V
}_
i
}
\frac
{
\partial
\Delta
V
_
i
}{
\partial
h
_
i
}
\left
(1 -
\frac
{
h
_
i
}{
n
_
d
\Delta
\tilde
{
V
}_
i
}
\frac
{
\partial
\Delta
\tilde
{
V
}_
i
}{
\partial
h
_
i
}
\right
)
^{
-1
}
.
\end{align}
In
\citet
{
hopkins2013
}
the lagrange multiplier is split into two parts such that
\begin{align}
\lambda
_
i =
\left
(
\frac
{
P
_
i
}{
\kappa
}
\frac
{
\Delta
\tilde
{
V
}^
2
_
i
}{
h
^{
n
_
d
}_
i
}
\right
)
\psi
_
i,
\label
{
eqn:sph:derivation:lambda
_
i
}
\end{align}
with
\begin{align}
\psi
_
i =
\frac
{
h
_
i
}{
n
_
d
\Delta
\tilde
{
V
}_
i
}
\frac
{
\partial
\Delta
V
_
i
}{
\partial
h
_
i
}
\left
(1 -
\frac
{
h
_
i
}{
n
_
d
\Delta
\tilde
{
V
}_
i
}
\frac
{
\partial
\Delta
\tilde
{
V
}_
i
}{
\partial
h
_
i
}
\right
)
^{
-1
}
.
\label
{
eqn:sph:derivation:psi
_
I
}
\end{align}
\subsection
{
Generalised Equation of Motion
}
Now that the lagrange multipliers have been constrained, the second half of the
equations in
\ref
{
eqn:sph:derivation:lmsum
}
(
$
q
_
i
=
\mathbf
{
r
}_
i
$
) can be used
to find the equation of motion. The differentials are given as
\begin{align}
\frac
{
\partial
L
}{
\partial
\dot
{
\mathbf
{
r
}}_
i
}
= m
_
i
\dot
{
\mathbf
{
r
}}_
i,
\quad
\frac
{
\partial
L
}{
\partial
\mathbf
{
r
}_
i
}
=
\sum
_{
j=1
}^
N P
_
j
\frac
{
\partial
\Delta
V
_
j
}{
\partial
\mathbf
{
r
}_
i
}
.
\end{align}
The differential of
\ref
{
eqn:sph:derivation:constraint
}
is
\begin{align}
\frac
{
\partial
\phi
_
j
}{
\partial
\mathbf
{
r
}_
i
}
=
-
\kappa
h
_
j
^{
n
_
d
}
\frac
{
\Delta
\tilde
{
V
}}{
\partial
\mathbf
{
r
}_
i
}
\frac
{
1
}{
\Delta
\tilde
{
V
}^
2
}
,
\end{align}
as
$
h
_
j
$
is an
\emph
{
independent coordinate
}
of a particle from
$
\mathbf
{
r
}_
j
$
implying that the partial differential
$
\partial
h
_
j
/
\partial
\mathbf
{
r
}_
i
=
0
$
.
We can substitute these into the equation of motion to find
\begin{align}
m
_
i
\frac
{
\mathrm
{
d
}
\mathbf
{
v
}_
i
}{
\mathrm
{
d
}
t
}
=
\sum
^
N
_{
j=1
}
P
_
j
\nabla
_
i
\Delta
V
_
j +
\lambda
_
j
\left
(
\kappa
\frac
{
h
_
j
^{
n
_
d
}}{
\Delta
\tilde
{
V
}_
j
^
2
}
\right
)
\left
(
-
\frac
{
\partial
\Delta
\tilde
{
V
}_
j
}{
\partial
\mathbf
{
r
}_
i
}
\right
).
\end{align}
This then gives the result from the substitution of the lagrange multipliers,
\begin{align}
m
_
i
\frac
{
\mathrm
{
d
}
\mathbf
{
v
}_
i
}{
\mathrm
{
d
}
t
}
=
\sum
^
N
_{
j=1
}
P
_
j
\nabla
_
i
\Delta
\tilde
{
V
}_
j
\psi
_
j +
P
_
j
\nabla
_
i
\Delta
V
_
j.
\label
{
eqn:sph:derivation:eom
}
\end{align}
Now we need to
\emph
{
specify
}
what we mean by volumes in terms of SPH
quantities. An example, familiar choice would be
$
\Delta
V
_
i
=
m
_
i
/
\rho
_
i
$
.
Notice that in this definition of volume there is one
\emph
{
particle-carried
}
property and one `field' or
\emph
{
smoothed
}
property. This turns out to be the
case for any smoothed property
\begin{align}
y
_
i =
\sum
^
N
_{
j=1
}
x
_
j W
_{
ij
}
(h
_
i)
\label
{
eqn:sph:derivation:smoothed
}
\end{align}
for a given particle-carried scalar
$
x
_
i
$
. This is also true for the
$
\Delta
\tilde
{
V
}_
i
=
\tilde
{
x
}_
i
/
\tilde
{
y
}_
i
$
. Using these specifications, we can write
down the volume differentials,
\begin{align}
\frac
{
\partial
\Delta
V
_
i
}{
\partial
h
_
i
}
=
-
\frac
{
x
_
i
}{
y
_
i
^
2
}
\frac
{
\partial
y
_
i
}{
\partial
h
_
i
}
,
\quad
\nabla
_
i
\Delta
V
_
j = -
\frac
{
x
_
i
}{
y
_
i
^
2
}
\nabla
_
i y
_
j.
\label
{
eqn:sph:derivation:volumediffs
}
\end{align}
The spatial differential is fairly straightforward and is given by
\begin{align}
\nabla
_
i y
_
j =
\nabla
_
i W
_{
ij
}
(h
_
j)
+
\delta
_{
ij
}
\sum
_{
k=1
}^
N
\nabla
_
i W
_{
ik
}
(h
_
i).
\label
{
eqn:sph:derivation:nablay
}
\end{align}
The differential with respect to the smoothing length is also straightforward
after remembering that
$
W
_{
ij
}
(
h
_
j
)
=
w
(
|r
_{
ij
}
|
/
h
_
j
)/
h
_
j
^{
n
_
d
}$
. Then,
\begin{align}
\frac
{
\partial
y
_
i
}{
\partial
h
_
i
}
= -
\sum
_{
j=1
}^
N
\frac
{
x
_
j
}{
h
_
j
}
\left
[
n
_
d W
_{
ij
}
(h
_
i) +
\frac
{
|r
_{
ij
}
|
}{
h
_
i
}
\left
.
\frac
{
\partial
W(u)
}{
\partial
u
}
\right
|
_{
u=
\frac
{
|r
_{
ij
}
|
}{
h
_
i
}}
\right
].
\label
{
eqn:sph:derivation:dydh
}
\end{align}
Finally, putting all of the above together, we can reach a
formulation-independent equation of motion for SPH,
\begin{align}
\frac
{
\mathrm
{
d
}
\vec
{
v
}_
i
}{
\mathrm
{
d
}
t
}
= -
\sum
_{
j=1
}^
N x
_
i x
_
j
&
\left
[
\frac
{
f
_{
ij
}
P
_
i
}{
y
_
i
^
2
}
\nabla
_
i W
_{
ij
}
(h
_
i)
\right
.
\nonumber
\\
&
\left
. + ~
\frac
{
f
_{
ji
}
P
_
j
}{
y
_
j
^
2
}
\nabla
_
i W
_{
ji
}
(h
_
j)
\right
],
\label
{
eqn:sph:derivation:spheom
}
\end{align}
with
\begin{align}
f
_{
ij
}
\equiv
1 -
\frac
{
\tilde
{
x
}_
j
}{
x
_
j
}
\left
(
\frac
{
h
_
i
}{
n
_
d
\tilde
{
y
}_
i
}
\frac
{
\partial
y
_
i
}{
\partial
h
_
i
}
\right
)
\left
(
1+
\frac
{
h
_
i
}{
n
_
d
\tilde
{
y
}_
i
}
\frac
{
\partial
\tilde
{
y
}_
i
}{
\partial
h
_
i
}
\right
)
^{
-1
}
\label
{
eqn:sph:derivation:fij
}
\end{align}
being the so-called `h-terms' that are essentially correction factors due to the
fact that
$
h
$
is allowed to vary.
Loading