Zeta function regularization (like any regularization procedure) is a choice about mathematics, since the limit can't depend on the procedure; the physics comes after the fact.
but how does an analytic continuation of the zeta function lead to a result that is compatible with experiments like the Casimir effect, do scientists know how this happened or not?
For the Casimir effect and related problems, yes, of course. The calculation is the difference between two infinite quantities, but the important point is that, while both quantities are infinite, their difference is finite. So it can be calculated by imposing any sort of regularization; zeta-function regularization is just one of many.
We do know, because we can show why the difference must be finite: It's because to define the infinite quantities it's necessary to introduce cutoffs, perform in this way the calculations in terms of finite quantities, and show that the results don't depend in how the cutoffs are defined.
The short answer is that the(ill-defined, since divergent) Σ w_n is defined as the value of the (well-defined, since convergent) Σ w_n f(w_n), where f is chosen so that the series converges. The non-trivial part of the process is the proof that the finite value is independent of the choice of the function f.
Changing s to z in the infinite sum, which gives the Riemann z-function, does not give the analytic continuity, automatically. For s=1, it looks smooth, at the plots. You can use s=1, perhaps, directly? I did some analysis, e.g. rewriting the sum into a halfinteger sum denoted sum(3/2), and then it factorise like
sum= (1-exp(-zln2))^-1(sum(3/2)+1)
so 'the thing' at s=0 is seen to appear in the denominator. Close to s=0, it fluctuates 'much', in numerical calculations, if you do not impose something else.