This question is really too broad. But there are some cool references that you can check. They are at the end of the answer.
The usual approach to represent realistic 3D fluids is to use a cell-based representation of the tracked portion of the world in which the fluid will be confined. Each cell has the velocity and pressure info, as shown in the image bellow. Instead of tracking the fluid in the world (Lagrangian viewpoint), the grid give us the possibility of look at fixed points in the space, the cells themselves (Eulerian viewpoint), which turns the problem easier to tackle.
Grid cell. p is the pressure, u, v and w are the velocity components for each axis.
To update the grid in each time step, the Navier-Stokes equation is used. The classic equation is a differential equation that represents the fluid advection, viscosity, pressure and external forces. For computational purposes, the equation is simplified so it does not have the viscosity factor (the error of the numerical solution can be seen as a viscosity component, however). This version of the equation is given bellow.
Momentum equation:
Incompressibility condition:
About the momentum equation, the left side and the first factor of the right side express the advection contribution, the second factor of the right side is the pressure contribution and the last factor of the right side is the external forces contribution.
Starting from the easier factor to compute, the external force contribution is just the acceleration caused by the net force of all fluid external forces (such as gravity, for example).
The advection contribution can be computed using the Semi-Lagrangian algorithm, which just uses the current velocity in a cell to travel to the past and get the value from another cell.
The pressure contribution is the most difficult, since usually it is also used to enforce the incompressibility condition and boundary conditions. It can be shown that the pressure solving leads to a Poisson equation, resulting in a linear system. To solve this equation the Modified Incomplete Cholesky Conjugate Gradient Algorithm can be used. It is fast but difficult to parallelise and implement. On the other hand, the Jacobi iterations are easy to implement and parallelise, so they are being used in shader implementations of fluid dynamics.
With this kind of representation the collision response for fluids comes practically for free, because of the advection and pressure updates. You just need to make sure that all fluid is inside the controlled world space.
In these references you can check shader implementations:
2D: GPU Gems 2D fluid.
3D: GPU Gems 3D fluid.
This last reference is a really great book that explains a lot about fluid simulation in detail: Robert Bridson's Fluid Simulation for Computer Graphics.