notes of SPH

Introduction

After no progress has been made in thermal-chemical nonequilibrium for months, I begin to doubt whether my ability can meet the need of basic scientific research work. So I put aside my current job, and spend a day reading a book[^book] about SPH method. Then I try the code provided by the author and gain much fun and confidence.
Here are some method and results.

[^book]: Smoothed Particle Hydrodynamics, a meshfree particle method G.R. Liu, M.B. Liu 2003

Method

The pith and marrow of SPH are fully embodied in that a function and its derivative can be represented in an integral form using dirac function or nascent delta function. It’s called kernel approximation in SPH.
$$
\begin{equation}
\left\langle f (\bf{x}) \right\rangle = \int\limits_\Omega f ({\bf{x’}}) W( {\bf{x}} - {\bf{x’}},h )d{\bf{x’}}
\end{equation}
$$

Then particle approximation is performed to conver the integral representation to the form of discretized particle approximation.

$$
\begin{equation}
\left\langle f\left( {\bf{x}}_i \right) \right\rangle = \sum\limits_{j = 1}^N \frac{m_j}{\rho_j} f\left( {\bf{x}}_j \right) \cdot W_{ij}
\end{equation}
$$

ODEs are produced when kernel approximation and particle approximation are applied to Navier-Stokes equations.

$$
\begin{equation}
\frac{d\rho_i}{dt} = \rho_i \sum\limits_{j = 1}^N \frac{m_j}{\rho_j}v_{ij}^\beta \cdot \frac{\partial W_{ij}}{\partial x_i^\beta }
\end{equation}
$$

$$
\begin{equation}
\frac{dv_i^\alpha }{dt} = \sum\limits_{j = 1}^N {m_j\frac{\sigma_i^{\alpha \beta } + \sigma_j^{\alpha \beta }}{\rho_i \rho_j}} \frac{\partial W_{ij}}{\partial x_i^\beta }
\end{equation}
$$

$$
\begin{equation}
\frac{de_i}{dt} =\frac{1}{2}\sum\limits_{j = 1}^N m_j\frac{p_i + p_j}{\rho_i \rho_j} v_{ij}^\beta \cdot \frac{\partial W_{ij}}{\partial x_i^\beta }+\frac{\mu_i}{2\rho_i}\varepsilon_i^{\alpha \beta} \varepsilon_i^{\alpha \beta}
\end{equation}
$$

Last time integration algorithm is achieved.

Result

I use the code provided by the book[^1] to simulate the shear cavity flow. Here is the result.

velocity distributions for the shear driven cavity flow

Some techniques

Fortran syntax

Write multi files for tecplot with time.

1
2
3
4
5
write(filename,'(i8)' ) itimestep
open (4, file = './data/' //trim(adjustl(filename))//'.dat')
write(4,*) "TITLE='CavityFlow'"
write(4,*) "VARIABLES = x,y,u,v"
write(4,*) 'ZONE AUXDATA time ="',time,'"'

Tecplot

Macro

References

1
2
3
4
5
6
$!VarSet |NumLoop|=100
$!Loop |NumLoop|
$!VarSet |num|=(|Loop|*500)
$!READDATASET '"D:\data\|num|.dat" '
$!EXPORTSETUP EXPORTFNAME = 'D:\data\png\|Loop|.png'
$!EndLoop

Dynamic text

1
&(auxzone[1]:time)

Gif movie gear