ACSE-4 Project 2: A Smoothed Particle Hydrodynamics Solver
Introduction to SPH
Smoothed Particle Hydrodynamics (SPH) is a meshless method for solving the Navier-Stokes equation. As it is meshless and Lagrangian (reference frame follows the fluid motion), it is ideal for solving problems involving fluid flow with interfaces and free surfaces. In this task you are going to write an SPH-based simulator and use it to solve a wave generation problem. The equations listed below include some which are purely there for background context and some which are required within the SPH simulation code. Those which you will need to implement are surrounded by a box.
The basic idea of SPH is that you model individual points which follow the fluid and are each associated with a given mass of the fluid, 𝑚𝑖. In order to carry out differentiation on them, the discrete values of any field (e.g., pressure) associated with the points are converted into a continuous function using a smoothing kernel (hence the smoothed in the name):
〈𝐴〉(𝐱)=∑𝑁 𝑚 𝐴𝑗 𝑊(𝐱−𝐱,h) 𝑗=1 𝑗𝜌𝑗 𝒋
Where 〈𝐴〉(𝐱) is the smoothed value of 𝐴 at location 𝐱, which is a vector, 𝐴𝑗 is the value of 𝐴 associated with particle 𝑗, 𝜌𝑗 is the density of particle 𝑗 and 𝑊(𝐱 − 𝐱𝒋, h) is a smoothing kernel which depends upon the vector 𝐱 − 𝐱𝒋 and a characteristic smoothing length h.
In SPH we are typically only interested in the values at the location of the points. We will thus use the following notation:
〈𝐴〉(𝐱)=〈𝐴〉=∑𝑁 𝑚𝐴𝑗𝑊(𝐫,h) 𝒊 𝑖 𝑗=1𝑗𝜌𝑗 𝒊𝒋
Where 𝐫𝒊𝒋 = 𝐱𝒊 − 𝐱𝒋. Also note that as this is a smoothing rather than an interpolation, in general 〈𝐴〉𝑖 ≠ 𝐴𝑖. At any point there is therefore a choice that can be made as to whether to use the smoothed or nominal value at that point.
The smoothing kernel must possess the property that it integrates to 1 when integrated over its area (or volume in 3D). The form of this smoothing kernel is important. It should have compact support (i.e. it should go to zero at a finite distance) so that only the neighbouring particles are involved in the calculations. It should also be smooth and, ideally, close to normal in shape. A number of different kernels have been proposed, but in this project I wish you to use the cubic spline:
Note that this kernel is symmetric and thus only depends on the distance between the points being considered. The prefactor in the equation will depend upon the number of dimensions and this equation is thus for 2D problems only.
If we want to use SPH to approximate the Navier-Stokes equation we need to have approximations for the partial derivatives of the variables of interest:
1 − 3 𝑞2 + 3 𝑞3
𝑖𝑓 0 ≤ 𝑞 ≤ 1
𝑖𝑓1≤𝑞≤2 𝑖𝑓 𝑞 > 2
10 2 4 𝑊(𝐫,h)=7𝜋h2{ 1(2−𝑞)3
Where 𝑞 = |𝐫| h
4
0
∇𝐴≈∑𝑁 𝑚𝐴𝑗∇𝑊(𝐫,h) 𝑖𝑗=1𝑗𝜌𝑗 𝒊𝒋
This is useful as ∇𝑊(𝐫𝒊𝒋, h) is known analytically, while all the other values are also known for each point. This means that calculating the gradient of a quantity is no more computationally expensive than calculating the average value at a point. For symmetric kernels such as we will be using in this project, we can further simplify this equation by using the following relationship:
∇𝑊(𝐫𝒊𝒋, h) = 𝑑𝑊 (|𝐫𝑖𝑗|, h) 𝐞𝒊𝒋 𝑑𝑟
Where 𝐞𝒊𝒋 = 𝐫𝑖𝑗 and 𝑑𝑊 (𝑟, h) is a known function of 𝑟 and h, which you will need to derive by |𝐫𝑖𝑗| 𝑑𝑟
differentiating the kernel (note that you must do this analytically, not numerically). Substituting this back into the previous equation gives:
∇𝐴≈∑𝑁 𝑚𝐴𝑗𝑑𝑊(|𝐫|,h)𝐞
𝑖 𝑗=1𝑗𝜌𝑗𝑑𝑟 𝑖𝑗 𝒊𝒋
The problem with this equation is that it is neither symmetric nor anti-symmetric (the contribution of particle 𝑖 on particle 𝑗 is not the same as the contribution of particle 𝑗 on particle 𝑖). (Anti)symmetry is desirable if we wish the system to explicitly conserve momentum and mass.
To achieve symmetric or anti-symmetric differentials we can use one of the following two relationships:
𝜌∇𝐴 = ∇(𝜌𝐴) − 𝐴∇𝜌
∇𝐴 = ∇ (𝐴) + 𝐴 ∇𝜌 𝜌 𝜌 𝜌2
Substituting the previous relationship into these equations results in the following approximations:
𝜌∇𝐴≈∑𝑁 𝑚(𝐴 −𝐴)𝑑𝑊(|𝐫 |,h)𝐞 𝑗=1 𝑗 𝑗 𝑖 𝑑𝑟 𝑖𝑗 𝒊𝒋
∇𝐴 ≈∑𝑁 𝑚 (𝐴𝑖 + 𝐴𝑗 ) 𝑑𝑊(|𝐫 |,h)𝐞 𝜌 𝑗=1𝑗𝜌𝑖2𝜌𝑗2𝑑𝑟𝑖𝑗 𝒊𝒋
Approximating the Navier-Stokes Equation using SPH
-Symmetric(notethe 𝐞 =−𝐞 ) 𝒊𝒋 𝒋𝒊
-Anti-symmetric
We now have all the tools required to approximate the Navier-Stokes equation and the associated continuity equations:
𝐷𝐯 = − ∇𝑃 + 𝜇 ∇2𝐯 + 𝐠 𝐷𝑡 𝜌 𝜌
𝐷𝜌 = −𝜌∇ ∙ 𝐯 𝐷𝑡
It should be noted that in a Lagrangian reference frame 𝐷 ≡ 𝜕 when the time derivatives are 𝐷𝑡 𝜕𝑡
carried out within the reference frame of differential volumes of the moving fluid (which in the SPH context is equivalent to the reference frame of the particles).
We have thus far not dealt with approximations for second derivatives, which are required to approximate the shear stress term. We could do this by taking the 2nd derivative of the kernel, but this is not a good approximation unless there are a large number of particles in the support region of the kernel, which is computationally expensive. An alternative is to use the following relationship:
∇2𝐯=∇∙∇𝐯≈∑𝑁 𝑚(1+1)𝑑𝑊(|𝐫|,h)∇𝐯∙𝐞
𝑖 𝑖𝑗=1𝑗𝜌𝑖2𝜌𝑗2𝑑𝑟𝑖𝑗 𝑖𝑗𝒊𝒋
We can make use of the fact that ∇𝐯𝑖𝑗 ∙ 𝐞𝒊𝒋 is the velocity gradient in the 𝑖𝑗 direction. This can be approximated as the difference in the velocity in the 𝑖𝑗 direction divided by the distance over which this occurs. Combining this approximation with the previous gradient approximations results in the following approximations for the Navier-Stokes and continuity equations:
Note that the sums in the above equations MUST NOT include the combinations where 𝑖 = 𝑗 as this will result in a division by zero! They also need only include neighbouring particles (those within 2h of the particle being considered. Note the dot product in the continuity equation.
These are the main equations that you will need to solve. You will notice that you currently have 3 equations (momentum equations in the 𝑥 and 𝑦 directions and the continuity equation) and 4 unknowns (velocity in the 𝑥 and 𝑦 directions, pressure and density) at each point. There is therefore still a missing equation required to close the system.
The way we will be closing the system is to assume that it is slightly compressible and to use an equation of state, which relates the pressure to the density. In this work we will use the Tait equation, which is a very stiff equation of state:
Where 𝐵 = 𝜌0 𝑐02. 𝜌0 is the initial or reference density of the fluid and 𝑐0 is a numerical speed of 𝛾
sound, which should NOT be the same as the actual speed of sound (the simulation would be prohibitively slow). The value of 𝑐0 should be chosen to be about an order of magnitude larger than the highest sustained speed in the system (for your work you can use 𝑐0 = 20 𝑚/𝑠, though you should test its impact on the results). To make sure that the equation is stiff enough use 𝛾 = 7.
Time stepping in SPH
The above equations allow a particle’s acceleration and rate of change of density to be calculated as a function of its neighbour’s positions, velocities and densities/pressures. We now need to be able to update the particle’s position, velocity and density, with the updated density being used to calculate a new pressure. To save complexity we will be doing this explicitly. We could use a simple explicit Euler scheme:
𝜕𝐯𝑁𝑃𝑃𝑑𝑊 𝑁11𝑑𝑊𝐯 𝒊≈𝐚𝑖=−∑𝑚𝑗( 𝑖 + 𝑗) (|𝐫𝑖𝑗|,h)𝐞𝒊𝒋+𝜇∑𝑚𝑗( + ) (|𝐫𝑖𝑗|,h) 𝑖𝑗
𝜕𝑡 𝑗=1 𝜌𝑖2 𝜌𝑗2 𝑑𝑟 +𝐠
𝜕𝜌 𝑁 𝑑𝑊
𝑖 ≈𝐷𝑖 =∑𝑚𝑗 (|𝐫𝑖𝑗|,h)𝐯𝑖𝑗 ∙ 𝐞𝒊𝒋
𝑗=1 𝜌𝑖2
𝜌𝑗2 𝑑𝑟
|𝐫𝑖𝑗|
𝜕𝑡
𝑗=1
𝑑𝑟
𝑃 = 𝐵 (( 𝜌 )𝛾 − 1) 𝜌0
𝐱𝐢𝑡+1 = 𝐱𝐢𝑡 + ∆𝑡 𝐯𝐢𝑡 𝐯𝐢𝑡+1 = 𝐯𝐢𝑡 + ∆𝑡 𝐚𝐢𝑡 𝜌𝑖𝑡+1 = 𝜌𝑖𝑡 + ∆𝑡 𝐷𝑖𝑡
This, though, is only first order accurate and is likely to be not very stable unless very small time steps are used. A scheme which is second order accurate in time is a predictor-corrector scheme:
Half-step
Full-step
𝐱𝐢𝑡+1/2 = 𝐱𝐢𝑡 + 0.5∆𝑡 𝐯𝐢𝑡
𝐯 𝑡+1 = 𝐯 𝑡 + 0.5∆𝑡 𝐚 𝑡 𝐢2𝐢𝐢
𝜌𝑖𝑡+1/2 = 𝜌𝑖𝑡 + 0.5∆𝑡 𝐷𝑖𝑡
𝐱̅ 𝑡 + 1 / 2 = 𝐱 𝑡 + 0 . 5 ∆ 𝑡 𝐯 𝑡 + 1 / 2 𝐢𝐢𝐢
𝑡+1 𝑡 𝑡+1 𝐯̅𝐢 2=𝐯𝐢 +0.5∆𝑡𝐚𝐢 2
𝜌̅𝑖𝑡+1/2 = 𝜌𝑖𝑡 + 0.5∆𝑡 𝐷𝑖𝑡+1/2
𝐱 𝑡 + 1 = 2 𝐱̅ 𝑡 + 1 / 2 − 𝐱 𝑡 𝐢𝐢𝐢
𝐯 𝐢 𝑡 + 1 = 2 𝐯̅ 𝐢 𝑡 + 1 / 2 − 𝐯 𝐢 𝑡 𝜌𝑖𝑡+1 = 2 𝜌̅𝑖𝑡+1/2−𝜌𝑖𝑡
The added complexity with the predictor-corrector scheme is that 2 values for each of the variables need to be stored (the initial value and the half-step value) and the neighbour searching needs to be based on the appropriate 𝐱 values.
We also need to choose a time step that will ensure stability. In this system the following time step constraints need to be satisfied:
∆𝑡𝐶𝐹𝐿=𝑚𝑖𝑛{h } |𝐯𝑖𝑗|
∆𝑡𝐹 = 𝑚𝑖𝑛 {√ h } |𝐚𝑖|
∆𝑡𝐴 = 𝑚𝑖𝑛 { h } 𝑐0√(𝜌⁄𝜌0)𝛾−1
∆𝑡 = 𝐶𝐶𝐹𝐿 𝑚𝑖𝑛{∆𝑡𝐶𝐹𝐿,∆𝑡𝐹,∆𝑡𝐴}
With an appropriate value of 𝐶𝐶𝐹𝐿 being in the range 0.1 – 0.3 (top end of the range is quicker, but less stable/accurate). Ideally a dynamic time step should be used, but as long as the speed of sound is high enough, it should dominate in the example you are required to solve. You can therefore use the following fixed time step for your initial testing:
Other required relationships
Another variable that you need to choose in the system is the relationship between the initial particle spacing, ∆𝑥 and the characteristic smoothing length, h. You can use:
∆𝑡 = 0.1 h 𝑐0
h = 1.3∆𝑥
Once you have the simulator working you can test the impact of changing this proportionality. The initial mass of your particles also depends upon the initial particle spacing:
Smoothing the Density/Pressure field
In SPH the integration of the continuity equation will result in not only variations in the density (and thus the pressure) based on the macroscopic velocity gradients, but also on local variations in particle spacing and velocity. This will eventually result in a very “rough” density and pressure distribution.
To get around this problem the density should be smoothed every 10 to 20 time steps. The following smoothing should be used:
Note that the sums MUST include the 𝑖 = 𝑗 combinations. Finding Neighbours
Virtually all of the above equations require summing over neighbouring particles. If we were to naively search over all possible neighbours for each of the points required in the calculation, then the overall complexity of the algorithm would be 𝑂(𝑁2) in the number of particles, which would be prohibitively expensive for all but the smallest problems.
We can get around this problem by exploiting the fact that the kernel has compact support and that we therefore only need to search for particles that are within 2h of the particle that we are considering. Considering that the particles will also be quite evenly spaced, the best way to do this search is using a linked cell method.
The way this works is that the entire domain is divided into a grid of 2h wide cells. A time step (or half time step in the predictor-corrector case) starts with each particle being assigned to their appropriate cell. This involves looping over all the particles and thus has a complexity of 𝑂(𝑁). The calculation for each particle involves looping over all the particles in its cell, together with all the particles in the neighbouring cells. Since the number of particles in these neighbouring cells does not depend upon the number of particles in the system, the calculations required to sum over all the neighbouring particles of a single particle is 𝑂(1), with the calculations for all the particles thus being 𝑂(𝑁). The overall algorithm is thus 𝑂(𝑁).
𝑚𝑖 = ∆𝑥2𝜌0
∑𝑛 𝑊(𝐫 ,h) 𝜌𝑖 = 𝑗=1 𝒊𝒋
∑𝑛 𝑊(𝐫𝒊𝒋,h) 𝑗=1 𝜌𝑗
Figure 1: Searching using linked cell method
In order to aid the development of your program you have been given an implementation of the linked cell search method. This implementation is comparatively basic and does not include the ability to use different half-step locations. It will, though, show you the basics of how to implement this method and form a basis on which you can build and improve.
Boundaries and initial conditions in SPH
Being a Lagrangian method, domain boundaries can be problematic to implement in SPH. An easy way to do this is to put a layer of fixed particles 2h wide around the edge of the whole domain. These particles should interact with other particles in the same way as the fluid particles and should have a density (and thus also pressure) that is allowed to evolve in the same way as that of other particles. Their position should remain fixed, though, and their velocity should remain zero.
Despite the boundary particles possessing a pressure, fluid particles may occasionally leak out of the domain. The easiest thing to do is to delete them, but this results in a loss of fluid mass. You can push them back, but you need to watch out for what the appropriate density to give them is (if you don’t they will clump as the local density will be higher than the average density assigned to the particle). The final alternative is to give the walls a repulsive force that gets high enough at short range that the particles never leak. This can be effective, but needs tuning and can result in instabilities if done incorrectly. Note that in order to ascertain if a particle has leaked you do need to know the location of the boundaries. Think about storing boundaries in their own data structures, as this will be especially useful for implementing arbitrary shaped boundaries.
The figure below gives the initial configuration that should be used for these simulations. In these simulations the following initial particle spacing is used:
∆𝑥 = 0.2 𝑚
Remember to make ∆𝑥 a variable that can be changed in the simulation configuration so that you can easily carry out a convergence test. A warning when placing the points is that the code is likely to crash if there are 2 points that are either on top of one another or close to being on top of one another. I would therefore recommend writing a function which checks for concurrent or near concurrent points and run this function after the points have initially been place, but before the simulation runs. Make sure that you only remove one of any two near concurrent points! Also remember that if a fluid particle is near concurrent to a wall particle it is the fluid rather than the wall particle that must be removed.
The simulations will also be of water and so the appropriate physical properties should be used (𝜌 =1000𝑘𝑔 and𝜇=0.001𝑃𝑎𝑠).
0 𝑚3
3m
3m
20 m
2m
10 m
Figure 2: Initial configuration for the simulation. Red are mobile fluid particles and blue are fixed boundary particles.
Code Structure
With your main simulation code you should allow the geometry, initial fluid layout and physical parameters to be entered in a relatively user friendly manner (this can be either via an input file or hard-coded). The simulation code should output the current state of the simulation at specified output time intervals (note that the output time interval will typically be much larger than the simulation time step). The simulation code must be written in C++.
There should be a separate program which loads the output files and performs the required visualisation and post-processing on them. This can be written in Python or C++. A useful way of visualising data is using Paraview. We have given you a short piece of code that is able to generate Paraview compatible files, though you are free to choose how you wish to visualise the outputs.
What is required of each team?
Tasks
1) Create an SPH simulator that can run based on the initial problem geometry and conditions described above. These simulations should run until 𝑡 = 30 𝑠 Start by implementing a simple forward Euler timestep method, for which you may need to use a small timestep.
2) Output particle data to file in such a form that it can be used for subsequent post- processing. Remember that the time steps will be quite small and so only produce outputs at a user set simulation time output interval.
3) Create a post-processing code to plot the particle data stored in the output files produced by the SPH simulator. The output from this code should be an animation of the simulation results, saved in a sensible format.
4) Improve your time-stepping method by implementing the predictor-corrector scheme described above. Note that you should save your functioning Forward Euler version as a separate file; you are not expected to have both time-stepping methods implemented in the same code.
5) Measure the wave height and location of the crest as a function of time. What is the sloshing interval across the domain and how does this compare to the expected shallow water wave
speed of √𝑔h0 where 𝑔 is gravity and h0 is the average initial water depth? Why might
there be a discrepancy?
6) Simulate different shaped domains, especially ones with non-vertical and horizontal walls.
Try and simulate a shoaling wave (e.g. running up onto a beach).
7) Carry out a convergence study in which the resolution (Δ𝑥) of the simulation is varied. You
will need to decide upon a suitable metric for the error (for example, timescale for several wave transits). What is the order of convergence?
Note that you are not expected to complete all of these tasks, though you should complete tasks 1-3 as a minimum. The SPH algorithm outlined in this description is a “bare bones” one and there are many ways in which the performance and stability can be improved. Extra credit will be given for enhancements to the code and so I would recommend looking in the literature for some of these improvements starting with the papers of J.J. Monaghan who invented the SPH method.
Recommendation for how to initially proceed
I would recommend that you quite rapidly decide upon a data structure for the particles. This can be based directly on the neighbour search code that you have been given for the Forward Euler time- stepping. Once the data structure has been agreed, you can split the team into one group to implement the numerical solver and another group to handle data output and rendering. Code to write the data to file and the post-processing code can be developed without a functioning SPH simulator, using just the initial positions of the particles, which can be defined by the code given. It is important to get both aspects developed at the same time as you can’t tell if the code is working correctly if you can’t see the results.
In Task 4 you will need to modify this data structure to include the need to store two values of position, velocity and density for each particle in order to implement the predictor-corrector time- stepping scheme, as well as modifying the neighbour search to use the appropriate particle positions. You do not need to modify the output routine or post-processing software for the predictor-corrector time-stepping as only the full timestep state needs to be outputted.
Technical requirements
You should use the assigned GitHub repository exclusively for your project Your main simulation software must be written in C++,
You are not allowed to import any external libraries for the main simulation code beyond the standard ones associated with C++ or the parallel libraries for MPI or openMP.
Your code should demonstrate its health and sustainability through a method such as automated testing. You may use the testing framework supplied in your starter repository or any alternative, but should use Travis as your CI engine.
Marking Scheme
Software (70 marks)
Functionality (30 marks)
This will be assessed based on the code’s ability to perform the tasks 1-6 outlined above. The required functionality therefore includes:
1) The ability to run successfully for 30 seconds of simulation time.
2) The ability to output the results to file.
3) A post-processor that can load the outputs and save them in a suitable format for rendering
as an video/animation
4) The implementation of a 2nd order accurate time-stepping scheme.
5) The extension of the post-processing to include the tracking of wave velocity.
6) The ability to simulate other domain geometries.
Performance (20 marks)
This mark will be assessed based on a combination of the correctness of the solution and the speed of the code. Included in the assessment of performance will be the convergence analysis carried out in task 7.
Sustainability (20 marks)
As with all software projects, you should employ all the elements of best practice in software development that you have learned so far. A GitHub repository will be created for your project to host your software. The quality and sustainability of your software and its documentation will be assessed based on your final repository and how it evolves during the week. Specific attention will be given to the following elements:
1) Installation and usage instructions
2) Documentation (via a suitable format).
3) Coding style
4) Quality and coverage of automatic testing framework
5) General repository usage
6) Licensing
Presentation (20 marks)
On Friday afternoon, you must demonstrate your software to staff during a 90 minute interactive session. Each assessor will ask you to demonstrate one of the following five elements of your software.
1) Describe the algorithm used to carry out the task, including any enhancements that you have made to the code. Particular focus should be given to how the particle containment was addressed, as well as stability/smoothness of the solution was ensured.
2) Demonstrate how the user would interact with the code, including setting the input conditions and post-processing the output.
3) Describe and discuss the results obtained, including how they compare with simplified analytic prediction for wave speed. Why might discrepancies exist?
4) Describe the quantification of the performance of the code by means of the convergence analysis. How does this convergence compare to the expected convergence? How else was the code performance quantified and assessed?
5) How did you address the problem of flow in arbitrary shaped domains? How were the boundaries specified? How was both the placement of the boundary particles and the containment of the fluid particles addressed?
Teamwork (peer assessment; 10 marks)
After the presentations, you will complete a self-evaluation of your group’s performance. This will
inform the teamwork component of your mark. Please refer to the ACSE-4 guidelines for more
information about the assessment of teamwork.