I've been thinking about divergent series on and off, so maybe I could chip in.
Consider a sequence of numbers (in an arbitrary field, e.g. real numbers) $\{a_n\}$. You may ask about the sum of terms of this sequence, i. e. $\sum a_n$. If the limit $\lim_{N\rightarrow\infty} \sum^N |a_n|$ exists then the series is absolutely convergent and you may talk about the sum $\sum a_n$. In case the limit does not exist but $\lim_{N\rightarrow\infty} \sum^N a_n$ exists then the sequence is conditionally convergent, and as (I assume) Carl Witthoft commented above there is a theorem stating that you may sum the sequence in a different order and get a different result for the limit. In fact by judiciously rearranging you may get any number desired. I included this just to mention that although divergent series may seem most bizarre, in the sense of summing terms and that by each term it gets nearer a limit, only the absolutely convergent series make connection with our intuiton. So we may ask about making sense of series in general.
As G. H. Hardy's "Divergent Series" explain in page 6, the trick is to understand that our usual notion of sum of a serie is a way to define something we call a "sum". In other words given a sequence we have a map that atributes to this sequence a number. The "sum map" being the trivial operation of summing the terms if the series is absolutely convergent. The idea behind divergent series is to realize that this map althogh in a sense canonical is not unique.
To be more specific, consider the space $V$ of all sequences together with operations of addition and scalar multiplication (given two sequences $\{a_n\}$ and $\{b_n\}$ and a number $\lambda$ we define addition by $\{a_n\}+\{b_n\}=\{a_n+b_n\}$ and scalar multiplication by $\lambda\cdot\{a_n\}=\{\lambda a_n\}$). Now the space of sequences with these operations is a (infinite-dimensional) vector space (there is a good question about coordinates, since I am assuming one specific basis here, but let's not worry about this now). The absolutely convergent series can be seen to form a subspace $U$ of $V$ and the "sum" $S$ is just a linear functional on this subspace, $S:U\rightarrow\mathbb{R}$. The problem is that this functional is not defined anywhere else.
So to make sense of divergent series one asks if is there another map $S'$ defined on a subspace $W$ that contains $U\subset W$ such that when restricted to this subspace is just the usual sum, i.e. is there $S'$ such that $S'|_U=S$?
And in fact many such functionals do exist. And each of these we call a different "summation method" in the sense that it attributes a value to a sequence and that when such sequence corresponds to a convergent series it gives the usual values.
For instance, Cesaro Summation says that maybe the series does not converge because it keeps oscillating (like the $1-1+1\cdots$ you mentioned). Then we could take the arithmetic mean of the partial sums $s_n=a_1+\cdots a_n$ and define the Cesaro sum of the sequence $\{a_n\}$ as $\lim_{N\rightarrow\infty}\frac{1}{N}\sum^N s_n$. Is not hard to see that this gives the usual result for convergent series (although a bit obscure that it is a linear map), but it also gives $1/2$ to the alternating series of $1-1+1\cdots$. So one must give a new meaning to the word "sum", and then you can get new results. For instance Fejer's theorem roughly states that (given mild conditions) the Fourier series of a function may not be stricty convergent, but it is always summable in the sense of Cesaro. So it tells you that the worst divergence that appears in Fourier series is of oscillating type, i. e. the series never diverges to $\pm\infty$. Furthermore by Cesaro summing you can tell around which value the series oscillates about. But this does not "sum" the series in the sense of making it convergent in the usual sense.
Other ideas for functionals is by analytic continuation. The most obvious is the geometric series $\sum x^n=\frac{1}{1-x}$. Is is only convergent in the $|x|\leq 1$ radius, but one may use analytic continuation to turn around the problem at $x=1$ and say that $1+2+4+\cdots=\frac{1}{1-2}=-1$. In this case is easier to picture the linear algebra idea. The space $U$ is of all geometric series with $|x|\leq 1$ and $S$ is the usual sum. Now we introduce a functional $S'$ by $S':W\rightarrow\mathbb{R},x\mapsto\frac{1}{1-x}$ so that $W$ is now every geometric series except the one with $x=1$. A functional which reduces to this case for the geometric series but does it also for other power series is Abelian Summation.
So the idea is not to "sum" the series really (only absolutely convergent series sum in usual sense) but to redefine the notion of sum by generalizing the concept and then using this different notion to attribute a finite value to the serie through the corresponding sequence. This finite value should tell you something about the series, like the Cesaro sum tells you around which value the series oscillates or like Abel sum is able to reconstruct the function that generated the series. So summation methods are able to extract information from divergent series, and this is how to make sense of them.
With respect to the physics, is important to stress that perturbative series in quantum field theory are (generally) divergent but neither renormalization nor regularization have to do (fundamentally) with "summing" divergent series (zeta regularization being one technique, not mandatory, although useful). Rather what does occur is that it sometimes one gets a asymptotic serie in perturbation theory. In general there are different functions with the same asymptotic serie, but with supplementary information it may, or may not, occur that one can uniquely find the function with that specific asymptotic series. In this case one can use a summation method known as Borel Summation to fully reconstruct the entire function. When such thing happens in QFT is normally associated with the presence of some sort of instanton. You can take a further look in S. Weinbergs "Quantum Field Theory Vol. 2", page 283. So the idea is to get non-perturbative information out of the perturbation series, and not to tame some sort of infinity. Renormalization is something completely different (and much worse since in fact it is highly non-linear for starters).
For further information try finding a copy of Hardy's book (it's a gem), or for the linear algebra babble J. Boos, F. P. Cass "Classical and Modern Methods In Summability".
A little bit of pedantry first: NaCl has a face-centered cubic lattice structure with 2 atoms in the basis (1 Na and 1 Cl), which looks like a simple cubic structure only if we ignore that there are two different kinds of atoms. In the following, $r$ will be the minimum separation between atoms, i.e. the shortest distance between a sodium and a chlorine atom. In the case of NaCl, the equilibrium distance of $r$ is exactly half the lattice constant $r_0 = a/2$, and would be the lattice constant of the simple cubic lattice, if the two types of atoms were indeed identical. Going off your question, I assume this is what you used.
Regarding your calculation, Ashcroft and Mermin point out the convention that threw you off in a footnote, in my edition it's in chapter 20, ionic crystals, page 403, footnote 11, which reads
It is customary to calculate the cohesive energy per ion pair, rather than by ion. If there are $N$ ions, then there will be $N/2$ ion pairs.
Since your cube of side length $r$ only contains a total of half a cation and half an anion, that's only half an ion pair. The correctly sized volume to go with the given convention for $u$ is therefore one that contains a full ion pair, which must then consequently be twice as large, i.e. $V=2r^3$.
Best Answer
This is not a rigorous answer to OP's questions, as we have not done any calculations or estimates.
The Madelung constant for the $NaCl$ crystal (of infinite size) is formally the sum
$$\tag{1} \sum_{\vec{r}\in \mathbb{Z}^3\backslash\{\vec{0}\}}\frac{(-1)^{x+y+z}}{||\vec{r}||_2} .$$
Here the $p$-norm is defined as
$$\tag{2} || \vec{r} ||_p ~:=~ \sqrt[p]{|x|^p+|y|^p+|z|^p}, \qquad || \vec{r} ||_{\infty}~:=~\max(|x|,|y|,|z|). $$
The procedure of arbitrarily truncating the summation (1) after completing a cube or ball of a certain size $a$ (possibly with an extra condition of only allowing electrically neutral sizes $a$), and then letting the size $a\to \infty$, is from a physical perspective rather crude.
Such truncation method seems ripe for systematic errors. If such truncation method happens to produce the correct answer, it seems more like an accident than a reliable scientific method.
Intuitively, it seems better to introduce a sufficiently smooth (or at least continuous) regulator function $\eta:[0,\infty[\to [0,\infty[$ with
$$\tag{3} \eta(0)~=~1\qquad \text{and}\qquad \lim_{r\to\infty}\eta(r)~=~0,$$
next consider the sum
$$\tag{4} \sum_{\vec{r}\in \mathbb{Z}^3\backslash\{\vec{0}\}}\frac{(-1)^{x+y+z}}{||\vec{r}||_2}\eta( ||\vec{r}||_p ~\varepsilon), $$
and finally removing the regulator $\varepsilon\searrow 0^{+}$.
In the spirit of Terence Tao blog entry, one would expect that the regularized sum (1) would not depend on the regulator function $\eta$ (or the $p\geq 1$ in the $p$-norm) once the regulator function $\eta$ belongs to a sufficiently nice class, and that such definition would agree with the (somewhat un-intuitive) definition by Borwein, Borwein and Taylor via analytic continuation.