Electronic – 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 – 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}}&&Z_3 = Z_{C_3}\mid\mid \infty=Z_{C_3}\\\\ V_Y = V_X\,\frac{1}{1+\frac{R_2}{Z_2}}&&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}}&&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&=0\\6\:\omega\:\tau&=\left(\omega\:\tau\right)^3\\6&=\left(\omega\:\tau\right)^2\\\omega&=\frac{\sqrt{6}}{\tau}\\\therefore\\f&=\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.