The packing of spheres (or non-spherical particles) within a geometry can be achieved quite reliably using the Discrete Element Method. I refer you to the following introduction article (written by me and co-workers) on the topic: https://onlinelibrary.wiley.com/doi/full/10.1002/cjce.23501
Briefly, the discrete element method (DEM) is a molecular dynamics inspired techniques where you calculate the motion of each individual particles by integration Newton's law of motion for each of them. The particle-particle and particle-wall contacts are handled by allowing for minute overlap between the particles (which are very small). These overlaps are then used within simple contact models. These models may be simple (linear spring) or more complex (hertzian spring + damping) and they have parameters which you have to tune to reproduce the mechanical behavior of your particles. These parameters are generally the static coefficient of friction or the rolling coefficient of friction.
DEM is relatively computationnaly intensive because you need to track all particles. However, using MPI, you can generally simulation a large amount of particles (over a million). There are also extensions of DEM to non-spherical particles.
There are some good open source codes such as :
DEM can even be coupled to CFD in what is referred to as CFD-DEM. You can look at the following journal article for a short introduction to it (written by my PhD student and I):
https://onlinelibrary.wiley.com/doi/full/10.1002/cjce.23686
You can look at the following online movies made by my group for an illustration. For pure DEM : https://www.youtube.com/watch?v=__5UGx4fQps
For CFD-DEM: https://www.youtube.com/watch?v=0IQo1Yq-_x8