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}
leftlangle f (bf{x}) rightrangle = intlimits_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}
leftlangle fleft( {bf{x}}_i right) rightrangle = sumlimits_{j = 1}^N frac{m_j}{rho_j} fleft( {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{drho_i}{dt} = rho_i sumlimits_{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} = sumlimits_{j = 1}^N {m_jfrac{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}sumlimits_{j = 1}^N m_jfrac{p_i + p_j}{rho_i rho_j} v_{ij}^beta cdot frac{partial W_{ij}}{partial x_i^beta }+frac{mu_i}{2rho_i}varepsilon_i^{alpha beta} varepsilon_i^{alpha beta}
end{equation}
$$

Last time integration algorithm is achieved.

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

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,'"'