- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Article Processing Charges ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Recently Accepted Articles ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Science and Technology of Nuclear Installations
Volume 2012 (2012), Article ID 351985, 9 pages
Development of a Secondary SCRAM System for Fast Reactors and ADS Systems
1Institute for Advanced Nuclear Systems, SCK.CEN, 200 Boeretang, 2400 Mol, Belgium
2Laboratory for Agricultural Machinery and Processing, Katholieke Universiteit Leuven, 30 Kasteelpark Arenberg, 3001 Heverlee, Belgium
3Institute of Mechanics, Materials and Civil Engineering, Catholic University of Leuven, 2 Place du Levant, 1348 Louvain-la-Neuve, Belgium
Received 15 November 2011; Accepted 14 February 2012
Academic Editor: Alberto Talamo
Copyright © 2012 Simon Vanmaercke et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
One important safety aspect of any reactor is the ability to shutdown the reactor. A shutdown in an ADS can be done by stopping the accelerator or by lowering the multiplication factor of the reactor and thus by inserting negative reactivity. In current designs of liquid-metal-cooled GEN IV and ADS reactors reactivity insertion is based on absorber rods. Although these rod-based systems are duplicated to provide redundancy, they all have a common failure mode as a consequence of their identical operating mechanism, possible causes being a largely deformed core or blockage of the rod guidance channel. In this paper an overview of existing solutions for a complementary shut down system is given and a new concept is proposed. A tube is divided into two sections by means of aluminum seal. In the upper region, above the active core, spherical neutron-absorbing boron carbide particles are placed. In case of overpower and loss of coolant transients, the seal will melt. The absorber balls are then no longer supported and fall down into the active core region inserting a large negative reactivity. This system, which is not rod based, is under investigation, and its feasibility is verified both by experiments and simulations.
One of the most important safety features of all reactor types is the ability of shutting down under all circumstances. This is in particular true for GEN IV reactors because they are designed to be safer than currently existing reactors thus also the ability to shut down the chain reaction must also be more reliable. In current PWR reactors shutting down the reactor can be accomplished in two completely independent and diverse methods. The first method is the insertion of safety and control bars by means of gravitation or other passive methods. The second method is the dissolution of the neutron absorbing boric acid into the primary water. For the liquid-metal- and gas-cooled GEN IV reactors (LFR, SFR, and GFR), the second method cannot be used because there are no liquid absorbers that can be dissolved in sufficient quantity in the liquid metal or gas, and even if such an absorbent would exist, cleaning the liquid metal after a SCRAM would be very expensive. In a liquid-metal-cooled ADS, there are in principle two different ways to shut down the reactor. First the accelerator can be turned off, leading to a safe shutdown of the subcritical core. Second the reactor power can be lowered by decreasing the multiplication factor of the core by inserting negative reactivity. This is the same problem as inserting reactivity in a critical reactor. (Note that this only decreases the power level and thus is not as effective as in a critical reactor.) The problem is therefore relevant for both ADS and critical reactors.
This paper focuses on inserting negative reactivity in a liquid-metal-cooled core, using a diverse operating principle compared to control/saefty rods. First it gives a short overview of the existing solutions to this problem. It finally presents a new concept that can be used in liquid-metal-cooled reactors and replaces the dissolution of a neutron absorber into the primary coolant. This concept is demonstrated by both simulations and experiments.
2. Overview of Existing Concepts
2.1. Hydraulically Suspended Absorber Balls
The hydraulically suspended absorber balls concept [1, 2] was designed by Rockwell International for use in sodium-cooled reactors. An illustration of this concept is shown in Figure 1. It uses spheres made of a neutron-absorbing material, tantalum, which are hydraulically suspended by the upward flow of the sodium coolant. In case this flow is interrupted, either by loss of coolant or by loss of flow, this upward force disappears and the tantalum spheres drop into the active core region. In this way they insert a negative reactivity shutting down the fission chain reaction. This system has the advantage of being self-actuating but can also be activated by a SCRAM system since it has a valve that can shut down the flow in case of SCRAM.
2.2. Liquid Absorber with Melt Seal
This concept, called LIM (lithium injection module) and illustrated in Figure 2, uses liquid lithium neutron poison in tubes. It was conceived for use in the Rapid and Rapid-L [3, 4] lithium-cooled reactor design concepts. These are self-controlling reactor concepts which do not require the intervention of an operator. An aluminum seal keeps a liquid poison (Li6) above the active core region during normal operation. In case of overpower, loss of coolant, and loss of flow accidents, the seal temperature increases, and at approximately 650°C it melts, thus releasing the Li6 into the active core region. This system is fully self-actuating, but it cannot be activated by a SCRAM system. The disadvantage of this system is that Li6 is not a very good neutron absorber, is chemically reactive, and is quite expensive.
2.3. ALMR Ultimate Shutdown System
The American ALMR , advanced liquid metal reactor, also plans to have a complementary shutdown system. This system also uses absorber spheres that are kept in a central channel above the core in normal operation by means of a seal. In case of SCRAM, the seal is ruptured by means of a mechanical device. With the seal ruptured, the balls are then free to flow in the active core region. This system is not self-actuating and requires a SCRAM signal.
2.4. Gas Expansion Modules
A self-actuating system for use in gas fast reactors (GFRs) was proposed, based on the gas expansion modules . In this system a boron carbide rod is submerged into a liquid metal. The level of the liquid metal, and thus boron carbide rod, is controlled by the liquid metal vapor pressure in the system. In an anticipated transients without scram, temperatures will increase resulting in an increased vapor pressure, decreasing the level of the liquid metal and thus inserting the neutron absorbing boron carbide rod into the active core region. Although this system is self-actuating and even self-resetting, it has an operational similarity to the control rod-based systems.
2.5. Articulated Absorber Rod
A complementary SCRAM system using an articulated absorber was implemented in the German SNR-300 [7–9] reactor. The system consists of three absorber elements linked in a chain. When the reactor is scrammed, the chain is pulled up through the core by means of an accelerator spring. Due to the flexibility of the chain, the assembly will conform to a duct with a bow that is several times larger than the maximum bow calculated for the guide tube . However, other blocking scenarios may arise, and the articulated absorber rod might in those cases be too similar to a conventional absorber rod.
3. Presented Concept
The concept that is presented in this paper aims at combining some of the strong points of the currently existing concepts. The presented concept is shown in Figure 3 . It consists of a tube, with the same diameter as a fuel pin, in the central position of every fuel assembly. The tube is divided into two sections by means of aluminum seal similar to the melting seal in the LIM concept. In the upper region, above the active core, spherical neutron-absorbing boron carbide particles are placed. In case of overpower and loss of coolant transients, the seal will melt. The absorber balls are then no longer supported and fall down into the active core region inserting a large negative reactivity.
A solid absorber was selected here because liquid-absorbing materials are rare and their negative reactivity effect is much lower compared to solid absorbers such as tantalum and boron carbide. A solid absorber has however the disadvantage that the flow of solid particles is not as well understood as the flow of liquids. Therefore it is clear that the absorber particle flow will be an important aspect of study for this concept.
The concept is also self-actuating just like the hydraulically suspended absorber balls and the liquid insertion with melting seal. It is however mechanically simpler than the hydraulically suspended balls which might contribute to the overall reliability. Unlike the hydraulically suspended absorber balls, it is not resettable and the reaction times might be longer.
Because this system is placed in every assembly, it has a distributed character similar to the dissolution of boric acid, unlike the primary SCRAM system using absorber rods that are localized to a few fixed positions or the ALMR ultimate shutdown system that is situated only in the central channel. Due to the flow-like behavior of the absorber balls, it is less prone to failure due to channel deformation/blocking than the absorber rods even though this issue has to be investigated.
Given the newly presented concept, a few questions require an answer before it can actually be used.(i)Can it provide sufficient negative reactivity to shut down a reactor?(ii)Will the particle flow block by itself?(iii)What is the optimal diameter of the spheres?(iv)Will the particle flow block when in contact with the molten seal?(v)Will the system be fast enough?
The following sections will address most of these questions, although further research is still needed.
4. Neutronic Effectiveness
Two important requirements of a SCRAM system are first to be able to insert sufficient negative reactivity to shut down the reactor under all circumstances and second to have a minimal impact on the normal operation of the core.
These two requirements can be verified using static MCNPX  calculations. This has been done on a critical variant of the MYRRHA core . MYRRHA is a lead-bismuth-cooled experimental accelerator-driven system.
In the reference subcritical core of MYRRHA, four fuel assemblies were added in order to obtain a critical core. In every assembly one fuel pin in the center was replaced with a SCRAM pin. These pins have been tested for their effectiveness with several materials. For every material the reactivity in normal operating case, where the absorber spheres are above the critical core, was compared to the scrammed case, where the absorber case was in the critical core region. This value yields the reactivity worth of the SCRAM system.
Also the normal case, where the spheres are above the critical core, was compared to a reference case, without the SCRAM system. This yields the reactivity influence of the system at normal operation.
The density for the sphere stacking is 61% of the theoretical density of the material that was used. This is the packing factor found in experiments when spheres are poured into a random stacking and conform to the value described in the literature . Experiments have also shown that this packing factor is very constant, with maximal deviation of about 0.5% points.
Table 1 shows both the impact of the secondary SCRAM system on the normal operation, as well as the SCRAM worth. Most of the effect on normal operation is caused by removing the fuel as can be seen with very ineffective absorbers, such as natural lithium. By far the most effective absorber material is enriched boron carbide. Another interesting candidate is however tantalum; although its neutron worth is not very high, it has the advantage of being heavier than lead bismuth and would therefore work reliably even when the SCRAM channel would be flooded. As a first option, 90% enriched boron carbide is retained. When channel flooding proves to be an important issue, tantalum can be considered as a valuable alternative.
5. Absorber Particle Flow Study
The presented concept relies heavily on the flow of the absorber spheres to shut down the reactor. Since particle flows are less well understood than liquid flows, verification of the reliability of this principle is therefore vital. This verification has been done both experimentally as well as with simulations. In a first instance, the particle flow is simulated alone without interaction of the molten aluminum seal. In a later stage the simulation of the seal will be done together with the particles giving a full simulation of the secondary SCRAM system. In this section only the results concerning the particle flow are discussed.
5.1. Simulation of Absorber Particles Flow Dynamics
The simulation starts with the particles in an initial position. The next step is to find all particles that are in contact with each other and compute the overlap they have. Based on this overlap and other parameters such as contact history, the normal and tangential forces are computed using a contact force model. All the normal and tangential forces working on a particle are then summed to one force and moment working in the center of gravity of the particle. Finally the forces are integrated to a new velocity and the velocity to a new particle position, and the scheme is repeated.
The choice of the selected normal and tangential force models is largely determinant in the accuracy of the simulation.
The normal contact force model is a Hertz-Kuwabara-Kono model [15, 16]: With the contact stiffness and the viscous normal damping parameter of particle or , the overlap between particle and is given by , the effective radius is computed as The contact stiffness is given by with , , respectively, the Young modulus and the Poisson ratio.
These models very accurately describe the forces between spheres in contact and require only physical constants to be known, being the Young modulus, viscous damping constant, and the static and dynamic Coulomb coefficients of friction and .
Due to the number of particles that amounts to several hundred thousands and the very small time steps required, the computational cost of the simulations is very high. To be able to complete the simulations in a reasonable time a, GPU- (graphics processing unit) based version of the DEMeter code was developed and used. This allowed for a speedup of 30 compared to the CPU-based version.
5.2. Experimental Investigation of Absorber Particle Flow Dynamics
The simulations done with the Demeter tool are compared to experiments. The experimental setup is shown in Figure 5 and consists of a 6 mm cylindrical glass tube filled with spherical “Ballotini” glass particles with a diameter of μm. The tube diameter is representative for the diameter of fuel rods used in the MYRRHA design.
The glass particles are kept in their initial position by a magnetic solenoid valve. The particle velocities are measured at three different heights (22 mm, 78 mm, 130 mm) from the bottom section of the tube. After the tube has been filled, the entire setup is completely enclosed in a vacuum chamber. This geometry is different from the real safety system, in the sense that the lower part of the tube is not present in the experiment. This has been done because it is much easier to do the experiments with the solenoid valve with this geometry. Additionally it is not expected that the lower part of the tube will have a significant influence on the particle flow. To measure the velocities of the particles, an optical flow method  has been selected. This has the advantage that the measurement does not influence the particle flow at all since it is nonintrusive. Optical flow measurement is a technique where images of the particle flow are taken at high frame rates. After capture, the frames are processed and sub-consequent frames are compared to each other to determine the magnitude and direction of the motion.
5.3. Results of the Particle Flow Analysis
For a tube filled to a height of 200 mm with glass spheres, the velocity evolution found both by experiments and simulations is shown in Figure 6. The experiment has been repeated 20 times in order to obtain an error estimation on the experiment. The velocities for the simulations are a cross-sectional average of the particles that would be visible by the experimental sensor. On average 100 particles are taken into account for this averaging.
Surprisingly the experiments are highly repeatable, yielding very similar results. This can be seen in the generally small error bars plotted on Figure 6. The standard deviation on the total time needed to empty the tube is less than 3% This feature is of course an advantage for a safety system. The change in error bar range is caused by the fact that the optical sensor sometimes is unable to track sufficient particles to yield an accurate velocity value. During periods with reduced particle density, this will increase the error bar drastically.
The particle velocity profile shown in Figure 6 is a simple acceleration of the particles as would be expected for particles in free fall mode. However, the particles are not in free fall mode; they do have interaction with the tube wall. This can be seen by the fact that the acceleration of the particles is not 9.81 m/s2 but closer to 6 m/s2. These interactions thus cause an energy loss. This energy loss however appears to be quite independent of the particle velocity, proven by the straight line in the velocity profile of the particles.
At the end of the graph, it can be observed that the velocity suddenly drops back to 0 m/s. This is simply due to the fact that the tube has been emptied and that there are no more particles to measure. The simulations predict the experimental data quite well, with a maximal difference in predicted time needed to empty the tube of about 5%.
5.4. Ideal Particle Size
The diameter of the absorber spheres is one of the few remaining free parameters that still can be optimized. Other parameters such as tube diameter and tube wall thickness, are determined by the necessity of not influencing the normal thermal hydraulics of the reactor. The location of the seal is fixed at the top of the active core region. This is a logical choice because at that location temperatures will be the highest, leading to shorter reaction times.
The absorber sphere size can be optimized for several goals: firstly it should be the size with a very low probability for blocking, next it should flow as fast as possible, and finally a high packing factor is desirable since it will increase the neutronic efficiency.
When the particle size becomes too small ≤100–250 μm, it is no longer the gravity force combined with collision forces that dominate and determine the movement of the particles. Other forces such as electrostatic, air friction, and especially Van der Waals forces become dominant. The latter one is an attractive force and reduces the free flowing capabilities of the particles. On top of deteriorated flow properties and a seriously increased risk of flow blockages, the flow becomes less easily predictable by simulations, due to the difficult to model electrostatic, and Van der Waals forces.
Choosing particles that are very large such as 1/2 or 1/3 of the tube diameter is also disadvantageous for several reasons. First of all they have a relatively high risk of blocking. Simulations and experiments show that the risk of blocking might be as high as 90% for diameters larger than 1.8 mm, although this risk seems to be very sensitive to the diameter in this region. Additionally the flow with these particles is slower compared to smaller sphere flows as can be seen in Figure 8 where the particle velocity is plotted for several diameters. Finally they have a lower packing factor as can be seen in Figure 7 showing the packing fraction as a function of the diameter. The packing factor is a decreasing function of the particle size, which can be expected since the wall effect on the stacking becomes increasingly important.
For the proposed complementary SCRAM system, this limits the possible range for the particle size between 500 μm and 1000 μm or between 1/12th and 1/6th of the tube diameter. The exact choice here is not very important; it can be seen in Figure 8 that in this diameter range the particle flow velocities are very insensitive to a change in particle diameter. Also the packing factor is not very sensitive to the exact choice of the diameter as can be seen in Figure 7 where the difference in packing factor between 500 μm and 1000 μm is less than 3%. For the envisioned range of particle diameters between 500 μm and 1000 μm, it was found impossible to find a single simulation that showed a blocking flow in over 3000 simulations. Also experimentally using 500 μm spheres, not a single experiment showed a blocking flow in over 400 experiments.
6. Absorber Particle-Metal Seal Interaction
This section describes the modeling of the interaction of the absorber particles and the aluminum seal that keeps the particles above the seal during normal operation and melts during a transient.
6.1. Molten Seal Modeling
The melting dynamics of the aluminum seal is modeled using smoothed particle hydrodynamics (SPH) . SPH is a Lagrangian particle method that does not require a computational mesh. It can be used to model a compressible fluid moving arbitrarily in three dimensions. SPH is very suitable to model free surface flows , and due to its particle-based nature, interaction with solid DEM particles can be easily implemented, which is the reason why this technique was selected to simulate the melting seal. The SPH equations for a liquid are obtained from the continuum equations of fluids dynamics by interpolating density, velocity, and so forth, from a set of points which may be disordered . The standard smoothed particle hydrodynamics method is extended with an energy transport equation , a simple melt model for the seal, and a surface tension model .
It was found experimentally that aluminum seals with a thickness of 1 mm and a diameter of 6 mm do not flow when molten but are kept steady by the surface tension. The seal shape deforms slightly as can be seen in Figure 9 where the experimentally obtained shape of the molten seal is compared to the predicted shape by the simulation. In order to have a functional secondary SCRAM system, the seal thickness has been reduced to 0.25 mm.
6.2. Seal-Particle Interaction
The aluminum seal has mechanical interaction with the solid absorber spheres both before the activation, when the aluminum seal is supposed to keep the particle above the core, and during the activation, when liquid aluminum seal can possibly mix with the solid particles. During the possible mixing between the solid particles and the molten seal, the seal could resolidify and form a blocking particle-seal conglomerate.
Because no wetting interaction between the aluminum and the boron carbide was observed, it was deemed sufficient to model the interaction using a Lennard-Jones potential barrier : where is a strength parameter, scaling the Lennard-Jones force, is the distance between points and , and is the distance at which the force is exactly zero, sometimes also called repulsion distance. If the distance between the two particles is shorter than , there is a repulsive force between the two particles; when the distance between the two particles is larger than , there is an attractive force that quickly drops off with increasing distance.
The combined particle/seal simulation with a 0.25 mm seal is shown in Figure 10. The outer surface of the secondary SCRAM system is supposed to have an instantaneous temperature increase from 450°C to 700°C or 50°C above the melting temperature of the aluminum seal. As can be seen, the particles rupture the liquid seal and are released into the core. The release of the particles, shown in the second image in Figure 10, occurs after 450 ms. The seal remains attached to the tube itself during the release of the particles and remains there afterwards.
In order to validate these results, the combined experiment has been repeated with a 0.25 mm seal. An infrared image sequence of the experimental particle release is shown in Figure 11. For better comparison, the simulation with the 0.25 mm seal, shown in Figure 10, is shown again in bottom view in Figure 12. These experiments have been repeated 5 times with consistent results.
The fact that the seal remains attached to the tube both in experiments and simulations reduces the interaction probability of the particles and the molten seal and therefore increases the overall reliability of the system.
In this paper the need for a secondary method to insert negative reactivity in an ADS or critical variant of an ADS prototype was explained. The currently proposed solutions for a reactivity insertion system have been summarized, and a new concept proposed. This concept uses spherical absorber particles that are placed in an empty tube with a similar diameter as a fuel pin. In normal operation they are kept above the active core region by means of an aluminum seal. In case of an accident, the temperature increases and melts the aluminum seal. This releases the particles into the active core region.
The flow of these particles is not as well understood as the flow of a liquid. Therefore this paper first addresses the flow of the absorber spheres, with both simulations as well as experiments. The flow was simulated with the discrete element method (DEM) and fits the experimental results well. Also, the ideal diameter range was determined to be between 500 μm and 1000 μm or between 1/12th and 1/6th of the tube diameter.
As a next step, the possible interaction between the solid absorber particle and the aluminum seal has been modeled. First the seal behavior has been simulated using the smoothed particle hydrodynamics (SPH) method, which has been extended with an energy transport equation, melt behavior, and surface tension model. This model has consequently been coupled to the particle simulation, and a combined simulation has been done. Using this model it was found that the time required to release the particles by melt of the seal in case that the outer surface of the system is 700°C is less than 500 ms. Also it was found that the liquid seal remains attached to the tube when molten and after release of the particles. This reduces the likelihood of interaction between the liquid aluminum and the particles and thus reduces blocking likelihood of the particles.
In the future unprotected and protected transients need to be done using a system code such as RELAP in order to verify that the proposed system is able to protect the reactor during these transients, without interfering with the primary SCRAM system.
- R. K. Paschall and A. S. Jackola, “Hydraulically supported absorber balls iss (inherent shutdown system)—water loop testing, absorber column,” Tech. Rep., Rockwell International, 1976.
- E. R. Specht, R. K. Paschall, M. Marquette, and A. Jackola, “Hydraulically supported absorber balls shutdown system for inherently safe lmfbr's,” in Proceedings of the International Meeting on Fast Reacctor Safety and Related Physics, vol. 3, p. 683, 1976.
- M. Kambe, H. Tsunoda, K. Mishima, and T. Iwamura, “RAPID-L operator-free fast reactor concept without any control rods,” Nuclear Technology, vol. 143, no. 1, pp. 11–21, 2003.
- M. Kambe and M. Uotani, “Design and development of fast breeder reactor passive reactivity control systems: LEM and LIM,” Nuclear Technology, vol. 122, no. 2, pp. 179–195, 1998.
- E. L. Gluekler, “U.S. Advanced liquid metal reactor (ALMR),” Progress in Nuclear Energy, vol. 31, no. 1-2, pp. 43–61, 1997.
- M. J. Driscoll, M. A. Pope, and P. Hejzlar, “Self-actuated reactivity insertion device for GFR service,” Transactions of the American Nuclear Society, vol. 89, pp. 573–576, 2003.
- OECD, “Status of LMFBR safety technology improving the performance and reliability of protection and shutdown systems,” OECD CSNI-R-69:0, 1983.
- K. Marten and H. Hoffmann, “Fluid dynamic investigations of snr-300 absorbers of the first and second shutdown unit,” Tech. Rep., Karlsruhe Nuclear Research Center, 1980.
- F. H. Morgenstern, F. Brandl, V. Ertel, and A. Schoensiegel, “The plant protection system of the SNR-300,” in Proceedings of the IEEE International Meeting on Fast Reactor Safety Technology, pp. 2602–2611, 1979.
- S. Vanmaercke, G. van den Eynde, E. Tijskens, and Y. Bartosiewicz, “Design of a complementary scram system for liquid metal cooled nuclear reactors,” Nuclear Engineering and Design, vol. 243, pp. 87–94, 2012.
- D. B. Pelowitz, MCNPX User Manual, RSICC, 2005.
- G. Van Den Eynde, E. Malambu, H. A. Abderrahim et al., “Neutronic design of the XT-ADS core,” in Proceedings of the International Conference on the Physics of Reactors (PHYSOR '08), vol. 3, pp. 2110–2115, 2008.
- F. A. Dullien, Porous Media: Fluid Transport and Pore Structure, Academic Press, 1991.
- P. Van Liedekerke, E. Tijskens, E. Dintwa, J. Anthonis, and H. Ramon, “A discrete element model for simulation of a spinning disc fertilizer spreader,” Powder Technology, vol. 170, pp. 348–360, 2009.
- N. V. Brilliantov, F. Spahn, J. M. Hertzsch, and T. Pöschel, “Model for collisions in granular gases,” Physical Review E, vol. 53, no. 5, pp. 5382–5392, 1996.
- G. Kuwabara and K. Kono, “Restitution coefficient in a collision between 2 spheres,” Japanese Journal of Applied Physics Part 1, vol. 26, no. 8, pp. 1230–1233, 1987.
- P. Franken, S. Francois, E. Tijskens, and G. Degrande, “A tangential force-displacement model for elastic frictional contact between the particles in triaxial test simulations,” in Proceedings of the International Conference on Particle-Based Methods Fundamentals and Applications, Barcelona, Spain, 2011.
- L. Vu-Quoc and X. Zhang, “Accurate and efficient tangential force-displacement model for elastic frictional contact in particle-flow simulations,” Mechanics of Materials, vol. 31, no. 4, pp. 235–269, 1999.
- R. D. Mindlin and H. Deresiewicz, “Elastic spheres in contact under varying oblique forces,” ASME Journal Applied Mechanics, vol. 20, pp. 327–344, 1953.
- B. D. Lucas and T. Kanade, “An iterative image registration technique with application to stereo vision,” in Proceedings of Imaging Understanding Workshop, pp. 121–130, 1981.
- V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot, “Meshless methods: a review and computer implementation aspects,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 763–813, 2008.
- M. Becker and M. Teschner, “Weakly compressible sph for free surface flows,” in Proceedings of the Eurographics/ACM SIGGRAPH Symposium on Computer Animation, 2007.
- J. J. Monaghan, “Smoothed particle hydrodynamics,” Annual Review of Astronomy and Astrophysics, vol. 30, no. 1, pp. 543–574, 1992.
- J. J. Monaghan, H. E. Huppert, and M. G. Worster, “Solidification using smoothed particle hydrodynamics,” Journal of Computational Physics, vol. 206, no. 2, pp. 684–705, 2005.
- A. Tartakovsky and P. Meakin, “Modeling of surface tension and contact angles with smoothed particle hydrodynamics,” Physical Review E, vol. 72, no. 2, Article ID 026301, pp. 1–9, 2005.
- J. E. Jones, “On the determination of molecular fields. from the equation of state of a gas,” Proceedings of the Royal Society of London Series A, vol. 106, pp. 463–477, 1924.