How to calculate the frequency of oscillation of this low-pass phase shift oscillator?
Your schematic uses a low pass filter approach. One way to mentally consider a solution is to look at the following schematic:



simulate this circuit &ndash; Schematic created using CircuitLab
Above, you can see that \$R_3\$ and \$Z_3\$ form a voltage divider that divides \$V_Y\$ into \$V_\text{OUT}\$. Also, \$R_2\$ and \$Z_2\$ form a voltage divider that divides \$V_X\$ into \$V_Y\$. Finally, \$R_1\$ and \$Z_1\$ form a voltage divider that divides \$V_\text{IN}\$ into \$V_X\$. It follows that:
$$\begin{align*}
\text{Each Stage}
\left\{
  \begin{array}{rl}
V_\text{OUT} = V_Y\,\frac{1}{1+\frac{R_3}{Z_3}}&amp;&amp;Z_3 = Z_{C_3}\mid\mid \infty=Z_{C_3}\\\\
V_Y = V_X\,\frac{1}{1+\frac{R_2}{Z_2}}&amp;&amp;Z_2 = Z_{C_2}\mid\mid \left(Z_3 + R_3\right)\\\\
V_X = V_\text{IN}\,\frac{1}{1+\frac{R_1}{Z_1}}&amp;&amp;Z_1 = Z_{C_1}\mid\mid \left(Z_2 + R_2\right)
  \end{array}
\right.
\end{align*}$$
$$\therefore \frac{V_\text{OUT}}{V_\text{IN}}=\frac{1}{1+\frac{R_1}{Z_1}}\cdot\frac{1}{1+\frac{R_2}{Z_2}}\cdot\frac{1}{1+\frac{R_3}{Z_3}}$$
Assuming \$R=R_1=R_2=R_3\$ and \$C=C_1=C_2=C_3\$ and setting \$\tau=R\:C\$ I get the following answer:
$$H\left(j\omega\right)=\frac{1}{1 - 5\:\left(\omega\:\tau\right)^2 + j\:\left[6\:\omega\:\tau-\left(\omega\:\tau\right)^3\right]}$$
For a phase shift of \$180^\circ\$, the imaginary part in the denominator goes to zero. So:
$$\begin{align*}6\:\omega\:\tau-\left(\omega\:\tau\right)^3&amp;=0\\6\:\omega\:\tau&amp;=\left(\omega\:\tau\right)^3\\6&amp;=\left(\omega\:\tau\right)^2\\\omega&amp;=\frac{\sqrt{6}}{\tau}\\\therefore\\f&amp;=\frac{\sqrt{6}}{2\pi\:R\:C}\end{align*}$$
Note that this differs from what you wrote.
I just tried out a simulation on LTspice using a decent R2R opamp (it can't handle more than about \$12.5\:\text{V}\$ between its rails.) Here's the results:

The period seems to be about \$5.8\:\text{s}\$. Which is close to the prediction.

Note
Per Tony's request, it's not difficult to work out the peak or peak-to-peak voltage for the resulting sine wave at the output.
Since the imaginary component of \$H\left(j\omega\right)\$ is 0, the magnitude is just \$\mid\:H\left(j\omega\right)\mid\:=\frac{1}{29}\$ (just plug in \$\omega=\frac{\sqrt{6}}{\tau}\$.) The RMS of the square wave at \$V_\text{IN}\$ is just \$6\:\text{V}\$ (as it is always either \$+6\:\text{V}\$ or else \$-6\:\text{V}\$.) So the output will be:
$$V_{\text{OUT}_\text{RMS}}=\frac{1}{29}\cdot \text{V}_{\text{IN}_\text{RMS}}=\frac{1}{29}\cdot 6\:{\text{V}_\text{RMS}}\approx 207\:{\text{mV}_\text{RMS}}$$
Since the output is a sine wave, the peak should be about \$\sqrt{2}\$ larger, or \$\approx 290\:\text{mV}_\text{P}\$.
As you can see, the image I included above shows an output peak that is close to this prediction, as well.