fluid dynamics - What is a good model for computing water dropping on a surface?

Before introducing my question I would like to outline the fact that I'm a coder, so I can be wrong when using some kind of terminology.

What is a good model for computing the flow, the shape and all the geometric properties of the water? ( and other similar liquids in general, like mixtures based on water or some more oily and viscous substances )

I assume that this kind of simulation will be based on particles, which means entities that are part of the domain that defines the "drop" of water, but what I need is the math, the equations and the mathematical model/s; I have no idea what is the math behind a liquid that meets a surface and if there is any algorithm or model that can be translated into machine code to compute that.

To keep things simple I'm not considering the air and what kind of properties the air can possibly modify, or any other extra variable outside the domain of this drop of liquid, like adding a inclined plane where the drops, let's say that for now I'm considering a vector and an X quantity of liquid, this quantity is than sprayed or just dropped on the surface and I should be able to reproduce a realistic dispersion of the water on the surface, even better I should be able to compute different liquids with different properties at the same time, like water and oil.

Answer

As I commented, I would think that any 3D hydrodynamics code would work. The basics of hydrodynamics can be summed up in the following five equations: \begin{eqnarray} \frac{\partial \rho}{\partial t}+\nabla\cdot\rho\mathbf{v}=0 \tag{1} \\ \frac{\partial \rho\mathbf{v}}{\partial t}+\nabla\cdot\left[\rho\mathbf{v}\otimes\mathbf{v}+P\mathbb I\right]=\rho\mathbf{g} \tag{2,3,4} \\ \frac{\partial E}{\partial t}+\nabla\cdot\left[\left(E+P\right)\mathbf{v}\right]=\rho\mathbf{g}\cdot\mathbf{v} \tag{5} \end{eqnarray} where $\rho$ is your mass density, $\mathbf{v}$ the velocity, $P=\left(\gamma-1\right)\left(E-\frac12\rho v^2\right)$ is the fluid pressure (assumes an ideal gas with adiabatic index $\gamma$), $\mathbb I$ the identity matrix (tensor),$E$ is the total energy (internal and kinetic), and $\mathbf{g}$ the gravitational acceleration. In order, these three describe mass conservation, momentum conservation and energy conservation.

Since computers are discrete objects, we define a volume $dx\cdot dy\cdot dz$ (usually $dx=dy=dz$ so the volume is $dx^3$, but it is not necessarily true) with a cell-centered position of $(x_i,\,y_j,\,z_k,\,t_n)$ where $i,j,k,n\in \mathbb Z$. Each variable is then defined by a volume-averaged value: $$\rho_{i,j,k}^n = \rho\left(x_i,\,y_j,\,z_k,\,t_n\right)$$

We can then numerically model the three conservation equations as (I'm going to use the mass conservation for the example) $$ \frac{\rho^{n+1}_{i,j,k}-\rho^n_{i,j,k}}{dt}=\frac{\pi_{x;\,i+1,j,k}^n-\pi_{x;\,i-1,j,k}^n}{2dx}+\frac{\pi_{y;\,i,j+1,k}^n-\pi_{y;\,i,j-1,k}^n}{2dy}+\frac{\pi_{z;\,i,j,k+1}^n-\pi_{z;\,i,j,k-1}^n}{2dz} $$ where $\pi=\rho\mathbf{v}$. The above uses a forward difference for the temporal derivative while the spatial derivatives use a central difference. There are some stability caveats with this equation (e.g. $dt\leq c_{s,max}dx/n_{dim}$ the time derivative is limited by the maximum wave speed divided by the number of dimensions times the cell length), but is fairly straight-forward to implement.

In three dimensions, you can either solve the equations directionally (i.e., $x,\,y,\,z$ solved independently but alternating order each time step) or you can combine them into what is called "corner transport" where the flux from say $\pi_{i-1,j-1,k}$ is accounted for in finding $\rho_{i,j,k}^{n+1}$. The latter choice is more difficult to code but provides a more accurate solution, while the former is pretty easy to implement.

For boundary conditions, you'd want a reflective boundary at the surface ($\mathbf{v}\cdot\hat{n}\to-\mathbf{v}\cdot\hat{n}$ where $\hat{n}$ is the normal direction to the surface) and probably extrapolated every where else ($\rho_{I,j,k}=\rho_{I-1,j,k}$ where $I$ is the maximum number of cells in the $x$ direction). These two boundaries and the above 5 equations should allow you to fully model the water droplet colliding with a surface.

You can also add viscous effects to Equations (2,3,4) by adding, to the RHS, the divergence of the viscosity tensor, $\tau$: $$ \tau=\rho\nu\left[\nabla\mathbf{v}+\left(\nabla\mathbf{v}\right)^T-\frac23\mathbb I\nabla\cdot\mathbf{v}\right] $$ with $\nu$ the kinematic viscosity. This complicates matters because you are introducing second-order derivatives ($\sim\nabla^2\mathbf v$) and the stability condition for these equations is $dt\leq \kappa dx^2/n_{dim}$ where $\kappa$ is the parabolic coefficient, in this case $\nu$. There are ways to get around this (e.g. using an implicit method that requires a matrix solver), but it definitely complicates matters.

---

With all of this, I conclude that you have two options:

1. write your own 3D code from scratch


2. search github for someone else's code that contains all the necessary physics and attribute the author when and where necessary

Given the difficulty in coding a multi-dimensional hydrodynamics code (personal experience here), it might be significantly easier for you to take Option 2, but my sole warning on that is this: the code you find and use is not a black box and should not be treated as such; you must understand what the code does and why before you can consider even running the code.