LAMMPS - Particle and fluid interaction - SPH

Hi all,

I am using LAMMPS for SPH. But my questions are general. I am trying to simulate a massive particle that would drop on a bed of water, so there can be 2 conditions: 1) the particle will displace some molecules of water resulting in splashes (as we see in the water collapse example) 2) the particle may diffuse with the water.

For the start, I am trying to simulate the first condition only. But my particle is just sliding over water molecules and water molecules do not displace at all.

My questions are:
i- I created water molecules in a region - is there any way I can remove the restriction of the region from water molecules so that they can move when hit by a particle?
ii- Can I specify the size/ diameter of the particle atoms?
iii- My last question is how can I make a read_data file to define the properties and initial position of the atoms/molecules.

I am sharing my input script here.



processors * * *

units		lj
atom_style	sph
boundary 	p p p

region 			box block 0 50 0 50 0 50 units lattice
create_box		2 box

## ------------- System definition ------------- ###

region			water cylinder z 20 20 10 INF INF

create_atoms	1 random 10000 	341341 water			#create 10000 atoms of type 1 in the region "water"
create_atoms	2 random 1 		127569 water			#create 150 atoms of type 2 in the region "mycylin"

### ------------- Simulation settings ------------- ###
mass			1 18										###attribute a mass of 18 (unitless) to atoms of type 1 
mass			2 60										###attribute a mass of 60 (unitless) to atoms of type 2 

group			water 		type 1
group			particle 	type 2

variable 			h 	equal 0.03
variable 			c 	equal 10.0 				# soundspeed for Tait's EOS
variable 			dt 	equal 0.1*${h}/${c} 	# upper limit of timestep based on CFL-like criterion

pair_style 			hybrid/overlay sph/taitwater/morris sph/rhosum 100 lj/cut ${h}
pair_coeff			* * sph/taitwater/morris 1000.0 ${c} 1.0e-3 ${h}  # use target density of 1000, soundspeed ${c} for Tait's EOS
pair_coeff 			* * sph/rhosum ${h}
pair_coeff 			* * lj/cut 0.005 ${h}

neigh_modify	every 5 delay 0 check yes

# time-integrate position, velocities, internal energy and density of inserted particles and water
fix             fix1 all sph

fix 			gfix particle gravity -9.81 vector 0 0 1 		# add gravity. This fix also computes potential energy of mass in gravity field.

#fix 			ID 		group-ID 	drag 	x 	y 	z 	fmag 	delta
fix 			dragfix particle drag 0 0 200 100 0 

compute            rho_peratom all sph/rho/atom
compute            e_peratom all sph/e/atom
compute            esph all reduce sum c_e_peratom
compute            ke all ke
variable           etot equal c_esph+c_ke+f_gfix

# adjust nevery | min. allowed dt | max. allowed dt | max. travel distance per dt # // ${dt} ~= CFL criterion 0.1*h/c
fix dtfix all dt/reset 1 NULL ${dt} 0.005 units box # use a variable timestep

dump 				dump_id1 all custom 100 combined.lammpstrj id type xs ys zs c_rho_peratom c_e_peratom fx fy fz
dump 				dump_id2 all xyz 100
dump_modify 		dump_id1 first yes
thermo 				100
thermo_style 		custom step ke c_esph v_etot f_gfix press time f_dtfix
thermo_modify 		norm no

timestep 			0.005
run 				10000

Many thanks

1 Like

There are no “molecules” in the SPH model, which is a grid-free continuum model.

Please see my comment above, what you can set instead are (local) density, energy and heat capacity.

Please see the documentation of the read_data command for a detailed description of the data file format and the number, order, and meaning of per-atom columns used atom style sph.

Please also make yourself more familiar with the theory governing SPH. Some general description is in Please check it out and the references in it. Please also note that this file was written for a much older version of LAMMPS (10+ years ago) and some names and technical details in LAMMPS and the SPH package have changed. The theory, however, has remained the same.

1 Like

Thanks for your response.
I have gone through the user guide, and yes some commands have changed that I noticed.
Still, I am trying to wrap my head around what is making my particles stationary even though I am integrating the whole range.
Is there any issue with the pair_style I am using?

Sorry, I don’t have the time to debug people’s inputs.

That’s alright. I completely understand.
Only started learning a few days back so I thought someone on the forum might help.


If you are new to this, you are moving ahead far to quickly. You should first study the example inputs line by line and learn what impact each of them has. Also, you should figure out which commands are required for the simulation and which are for producing output and analysis. Then you can move on to modifying them in small steps toward your actual simulation, again always checking that each modifications does what you expect.


Thanks a lot for the feedback.
Sure thing, I’m totally on it.
Hope to learn it soon.

Thanks again