Why is this series of square root of twos equal $\pi$?

It's the sequence of approximations obtained when you approximate the perimeter of the circle of diameter $1$ with inscribed regular $n$-gons for $n$ a power of $2$.

As I happen to have this TeXed' up, I'll offer:

Suppose regular $2^n$-gons are inscribed in a circle of radius $r$.

Suppose the side length(the length of one "face") $a_n$ of a the inscribed $2^{n}$-gon is known (so, $a_2$ is the side length of the square). To find the side length of the $2^{n+1}$-gon, one may apply the Pythagorean Theorem twice to obtain $$ \tag{1}a_{n+1} = r\sqrt{2-\sqrt{4-{a_n^2\over r^2}}} $$

Now, starting with a square, $$a_2=\sqrt 2 r.$$ Using the recursion formula (1) repeatedly gives: $$ a_3%= r\sqrt{2-\sqrt{4-{2r^2\over r^2}}} = r\sqrt{2-\sqrt2}, $$ $$ a_4%= r\sqrt{2-\sqrt{4-{ ( r\sqrt{2-\sqrt2})^2 \over r^2}}} = r\sqrt{2-\sqrt{4-{ ({2-\sqrt2} ) }}} = r\sqrt{2-\sqrt{{ {2+\sqrt2} }}}, $$ and $$ a_5%= r\sqrt{2-\sqrt{4-{ ( r\sqrt{2-\sqrt{{ {2+\sqrt2} }}} )^2\over r^2}}} = r\sqrt{2-\sqrt{ 2+\sqrt{{ {2+\sqrt2} }}} }. $$ $$\vdots$$

Let $b_n=2^n a_n$. Let $P_n=r\cdot b_n$ be the perimeter of the $2^n$-gon. Let $P$ be the perimeter of the circle. Then $$ \lim_{ n\rightarrow \infty} P_n = P. $$ Note that from the above identity, it follows that the ratio of the perimeter of a circle to its diameter must be a constant, namely $\lim\limits_{n \rightarrow \infty} b_n$. We call this number $\pi$.

---

Below are some particular calculations when the radius of the circle is $1/2$:

$$\eqalign{ P_2&=2^1\cdot\sqrt 2 \approx 2.82842712\cr P_3&=2^2\cdot\sqrt{2-\sqrt2}\approx 3.06146746\cr P_4&=2^3\cdot\sqrt{2-\sqrt{2+\sqrt2}}\approx3.12144515 \cr P_5&=2^4\cdot\sqrt{2-\sqrt{2+\sqrt{2+\sqrt2}}}\approx 3.13654849\cr P_6&=2^5\cdot\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt2}}}}\approx 3.14033116\cr P_7&=2^6\cdot\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt2}}}}}\approx 3.14127725\cr P_8&=2^7\cdot\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt2}}}}}}\approx 3.1415138 \cr P_9&=2^8\cdot\sqrt{2-\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt{2+\sqrt2}}}}}}}\approx 3.14157294 \cr } $$

---

For completedness:

Remark 1: Here is the proof that the recursion formula (1) holds:

Let $a_n$ be the side length of the $2^n$-gon.

To obtain the $2^{n+1}$-gon: take the "outer end point" of the radii of the circle that bisect the faces of the $2^n$-gon to form the new vertices of the $2^{n+1}$-gon.

We then have, for $a_{n+1}$, the scenario shown in the following diagram (not to scale):

Now $$ b^2=r^2-{a_n^2\over4}; $$ whence $$\eqalign{ a_{n+1}^2={a_n^2\over4} + \Biggl((r-\sqrt{ r^2-{a_n^2\over4}}\ \Biggr)^2 &={a_n^2\over4}+ r^2-2r\sqrt{r^2-{a_n^2\over4}}+r^2 -{a_n^2\over4}\cr &= 2r^2-2r\sqrt{r^2-{a_n^2\over4}}\cr &= 2r^2-r^2\sqrt{4-{a_n^2\over r^2}}\cr &= r^2 \Biggl(2-\sqrt{4-{a_n^2\over r^2}}\ \Biggr).}$$ And, thus $$ a_{n+1}= r \sqrt{2-\sqrt{4-{a_n^2\over r^2}}}. $$

Remark 2: To explain why limit $\lim\limits_{n\rightarrow\infty} P_n=P\ $ holds, I can do no better than refer you to Eric Naslund's comment in his answer.

See also, here.

Here is a slightly different way to see why it is the area of a $2^k$-gon. (It is really the same, I just also want to point out that the nested radical expression is the sin of $\frac{\pi}{2^k}$. This type of argument gives Vietas product formula for $\frac{2}{\pi}$)

Recall: If $a,b$ are sides of a triangle, and $\theta$ is the angle between them, then the area of this triangle is $\frac{1}{2}ab\sin(\theta)$.

For any regular $n$-gon inscribed in the unit circle, this means that by splitting it into the $n$ identical isosceles triangles with two sides equal to $1$, we have that $$\text{Area of regular n-gon}=\frac{n}{2}\cdot \sin\left(\frac{2\pi}{n}\right).$$

In the limit, this must approach $\pi$. Now, the nested radical expression can be rewritten as sin of angle as follows:

We have that $$\sin\left(\frac{\pi}{2^k}\right)=\frac{1}{2}\sqrt{2-\sqrt{2+\sqrt{2+\cdots+\sqrt{2}}}}.$$

Proof: Notice that $$\cos(x)=\sqrt{\frac{1}{2}\left(1+\cos(2x)\right)}.$$ Using this iteratively allows us to find $\cos\left(\frac{\pi}{2^k}\right)$. For example, since $\cos\left(\frac{\pi}{2}\right)=0$ we see that $$\cos\left(\frac{\pi}{4}\right)=\frac{1}{\sqrt{2}}=\frac{1}{2}\cdot \sqrt{2}$$ and similarly

$$\cos\left(\frac{\pi}{8}\right)=\sqrt{\frac{1}{2}\left(1+\frac{1}{\sqrt{2}}\right)}=\frac{1}{2}\cdot \sqrt{2+\sqrt{2}}.$$ Iterating again we have $$\cos\left(\frac{\pi}{16}\right)=\frac{1}{2}\cdot \sqrt{2+\sqrt{2+\sqrt{2}}},$$ and by induction we'll arrive at our result.

Since $\cos^2(x)+\sin^2(x)=1$, one final manipulation will allow us to conclude our result.

This is more of an extended comment than an answer, but I figure that since the OP is interested in the computational aspect ("...I wrote a little script to calculate it..."), this might be of some interest. The following is adapted from Ole Østerby's unpublished manuscript.

---

As already noted, the OP's limit is equivalent to saying that

$$\lim_{h\to 0}\frac{\sin\,\pi h}{h}=\pi$$

where we identify $h$ with $2^{-n}$. That is, to borrow notation from David's answer (and taking $r=1$):

$$b_n=2^n \sin\left(\frac{\pi}{2^n}\right)$$

If we expand $\dfrac{\sin\,\pi h}{h}$ as a series, like so:

$$\frac{\sin\,\pi h}{h}=\pi-\frac{\pi^3 h^2}{6}+\frac{\pi^5 h^4}{120}+\cdots$$

we see that only even powers of $h$ occur in the expansion (as expected, since the function in question is even).

If we halve $h$ (equivalently, increase $n$), we have a slightly more accurate approximation of $\pi$. The upshot is that one can take an appropriate linear combination of $\dfrac{\sin\,\pi h}{h}$ and $\dfrac{\sin(\pi h/2)}{h/2}$ to yield an even better approximation to $\pi$:

$$\frac13\left(\frac{4\sin(\pi h/2)}{h/2}-\frac{\sin\,\pi h}{h}\right)=\pi -\frac{\pi^5 h^4}{480}+\frac{\pi^7 h^6}{16128}+\cdots$$

This game can be repeatedly played, by taking successive linear combinations of values corresponding to $h/2$, $h/4$, $h/8$... The method is known as Richardson extrapolation.

(I'll note that I have already brought up Richardson extrapolation in a number of my previous answers, like this one or this one.)

More explicitly, taking $T_n^{(0)}=b_n$, and performing the recursion

$$T_j^{(n)}=T_{j}^{(n-1)}+\frac{T_{j}^{(n-1)}-T_{j-1}^{(n-1)}}{2^n-1}$$

the "diagonal" sequence $T_n^{(n)}$ is a sequence that converges faster to $\pi$ than the sequence $b_n$. Christiaan Huygens used this approach (way before even Richardson considered his extrapolation method) to refine the Archimedean estimates from circumscribing and inscribing polygons.

---

Sundry Mathematica code:

Table[2^(n - 1)*Sqrt[2 - Nest[Sqrt[2 + #1] & , 0, n - 2]] ==
      FunctionExpand[2^n*Sin[Pi/2^n]], {n, 2, 15}]

{True, True, True, True, True, True, True, True, True, True,
 True, True, True, True}

This verifies the equivalence of the iterated square root and sine expression.

Here is an implementation of the application of Richardson extrapolation to the computation of $\pi$:

huygensPi[n_Integer, prec_: MachinePrecision] := 
 Module[{f = 1, m, res, s, ta},
  res = {ta[1] = s = N[2, prec]};
  Do[
   If[k > 1, s = s/(2 + Sqrt[4 - s])];
   f *= 2; ta[k + 1] = f Sqrt[s];
   m = 1;
   Do[m *= 2;
    ta[j] = ta[j + 1] + (ta[j + 1] - ta[j])/(m - 1);, {j, k, 1, -1}];
   res = {res, ta[1]};, {k, n - 1}];
  Flatten[res]]

Note that I used a stabilized version of the recursion for generating the $b_n$ (f Sqrt[s] in the code) to minimize errors from subtractive cancellation.

Here's a sample run, where I generate 10 successive approximations to 25 significant digits:

huygensPi[10, 25] - Pi

{2.000000000000000000000000, 3.656854249492380195206755,
 3.173725640962868268333725, 3.13944246625722809089242,
 3.14157581875151427853903, 3.14159291451874033422144,
 3.14159265388327967181647, 3.14159265358866759077617,
 3.14159265358979303435864, 3.14159265358979323865872}

where the $10$-th approximant is good to 19 digits. By comparison, $b_{10}=2^{10}\sin\left(\dfrac{\pi}{2^{10}}\right)= 3.1415877\dots$ is only good to five digits.