Okay, after figuring out which paper you were trying to link to in the third link, I decided that it is better to just give an answer rather then a bunch of comments. So... there are several issues at large in your question. I hope I can address at least some of them.
The "big picture" problem you are implicitly getting at is the Hilbert problem of hydrodynamical limit of the Boltzmann equations: that intuitively the ensemble behaviour at the large, as model by a fluid as a vector field on a continuum, should be derivable from the individual behaviour of particles, as described by kinetic theory. Very loosely tied to this is the problem of global existence and regularity of Navier-Stokes.
If your goal is to solve the Navier-Stokes problem using the hydrodynamic limit, then you need to show that (a) there are globally unique classical solutions to the the Boltzmann equations and (b) that they converge in a suitably regular norm, in some rescaling limit, to a solution of Navier-Stokes. Neither step is anywhere close to being done.
As far as I know, there are no large data, globally unique, classical solutions to the Boltzmann equation. Period. If we drop some of the conditions, then yes: for small data (perturbation of Maxwellian), the recent work of Gressman and Strain (0912.0888) and Ukai et al (0912.1426) solve the problem for long-range interactions (so not all collision kernels are available). If you drop the criterion of global, there are quite a bit of old literature on local solutions, and if you drop the criterion of unique and classical, you have the DiPerna-Lions solutions (which also imposes an angular-cutoff condition that is not completely physical).
The work of Golse and Saint-Raymond that you linked to establishes the following: that the weak solution of DiPerna-Lions weakly converges to the well-known weak solutions of Leray for the Navier-Stokes problem. While this, in some sense, solve the problem of Hilbert, it is rather hopeless for a scheme trying to show global properties of Navier-Stokes: the class of Leray solutions are non-unique.
As I see it, to go down this route, you'd need to (i) prove an analogue of DiPerna-Lions, or to get around it completely differently, and arrive at global classical and unique solutions for Boltzmann. This is a difficult problem, but I was told that a lot of very good people are working on it. Then you'd need (ii) also to prove an analogue of Golse-Saint-Raymond in a stronger topology, or you can use Golse-Saint-Raymond to first obtain a weak-limit that is a Leray solution, and then show somehow that regularity is preserved under this limiting process. This second step is also rather formidable.
I hope this somewhat answers your question.
You are completely correct in your analysis of the structure one obtains by considering the
Schrodinger equation for Z electrons in a central potential due to a nucleus of charge Ze
when the Coulomb interaction between electrons as well as relativistic effects such as spin-orbit coupling are ignored. In such a model the periods would indeed be 2, 8, 18 etc. for
exactly the reasons you have described. In physics jargon the energy in this model depends only on the principal quantum number $n \in {\mathbb Z}, n>0$ and the allowed $\ell$ values
are $\ell\le n-1$. Orbitals with $\ell=0,1,2,3,4, \cdots$ are labelled by $s,p,d,f \cdots$ for
historical reasons. Thus at $n=1$ one can have one or two states in the $(1s)$ orbital (accounting for spin), at $n=2$ one has the $(2s)$ and $(2p)$ orbitals with $2(1+3)=8$
states. And at $n=3$ one should have $18$ states by filling the $(3s),(3p),(3d)$ orbitals.
But in the real world this is not what happens. This simple model does not give a correct description of atoms in the real world once you get past Argon. I believe the main effect leading to the breakdown is the Coulomb interaction between the electrons.
So there is no simple mathematical model based say just on the representation theory of $SU(2)$ and simple solutions to the Schrodinger equation which will account for the structure of the periodic table past Argon. However one could ask whether including Coulomb interactions between electrons does gives a model which correctly reproduces the next few rows of the periodic table past Argon. I am not an expert on this, but since I doubt there are physical chemists on MO I'll just give my rough sense of things.
To approach this problem with some level of rigor probably requires difficult numerical work
and my impression is that this is beyond the current state of the art. However there are
rough models which try to approximate what is going on by assuming that the interactions between electrons can be replaced by a spherically symmetric potential which is no longer
of the $1/r$ form. This leaves the shell structure as is, but can change the ordering of
which shells are filled first. In such a model instead of filling the $(3d)$ shell after
Argon one starts to fill the $(4s)$ and $(3d)$ shells in a somewhat complicated order. Eventually one fills the $(4s),(3d),(4p)$ shells and this leads to the line of the periodic
table starting at K and ending at Krypton.
Added note: There is one nice piece of mathematics associated with this problem that I should have mentioned, even if it doesn't by itself explain the detailed structure of the periodic table. When Coulomb interactions between electrons and relativistic effects are ignored the energy levels of the Schrodinger equation with a central $1/r$ potential depend only the quantum number $n$, but not on the quantum number $\ell$ which determines the representation $V_\ell$ of $SO(3)$ referred to above. When screening is included this is no longer the case and the energies depend on both $n$ and $\ell$. Why is this?
With a $1/r$ central potential there is an additional vector $\vec D$ which commutes with the Hamiltonian. Classically this vector is the Runge-Lenz vector and its conservation explains why the perihelion of elliptical orbits in a $1/r$ potential do not precess. Quantum mechanically the
commutation relations of the operators $\vec D$ along with the angular momentum operators $\vec L$ are those of the Lie algebra of $SO(4)$ (for bound states with negative energy). There are two Casimir invariants, one vanishes and the other is proportional to the energy.
As a result the energy spectrum depends only on $n$ and can be computed using group theory without ever solving the Schrodinger equation explicitly. Perturbations due to screening, that is from some averaged effect of the Coulomb interactions between electrons,
change the $1/r$ potential to some more general function of $r$ and break the symmetry
generated by the Runge-Lenz vector $\vec D$. As a result the energy levels depend on both $n$ and $\ell$.
Best Answer
To answer your first question:
Actually the assumption is not that the wave function and its derivative are continuous. That follows from the Schrödinger equation once you make the assumption that the probability amplitude $\langle \psi|\psi\rangle$ remains finite. That is the physical assumption. This is discussed in Chapter 1 of the first volume of Quantum mechanics by Cohen-Tannoudji, Diu and Laloe, for example. (Google books only has the second volume in English, it seems.)
More generally, you may have potentials which are distributional, in which case the wave function may still be continuous, but not even once-differentiable.
To answer your second question:
Once you deduce that the wave function is continuous, the equation itself tells you that the wave function cannot be twice differentiable, since the second derivative is given in terms of the potential, and this is not continuous.