The correct expression for the variance of the $i$th residual is explained here in detail.
I am using a slightly different notation but in the end it will all match with what you are working with.
Suppose we have a simple linear regression model $$y=\alpha+\beta x+\epsilon$$
where $\alpha+\beta x$ is the part of $y$ explained by $x$ and $\epsilon$ is the unexplained part or the error. Here $y$ is stochastic and $x$ is non-stochastic.
We consider the paired observations $(x_i,y_i)$ for $i=1,2,\ldots,n$ and assume that $\epsilon_i$ are i.i.d $\mathcal N(0,\sigma^2)$ for all $i$. This means we have some $Y_i\sim\mathcal N(\alpha+\beta x_i,\sigma^2)$, independently for all $i$.
Define $$s_{xx}=\sum_{i=1}^n (x_i-\bar x)^2\qquad,\qquad s_{yy}=\sum_{i=1}^n (y_i-\bar y)^2$$ and $$s_{xy}=\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)$$
From the normal equations, we have the least square estimates of $\alpha$ and $\beta$ :
$$\hat\alpha=\bar y-\hat\beta \bar x\qquad,\qquad\hat\beta=\frac{s_{xy}}{s_{xx}}$$
Let the residual variance be $$s^2=\frac{1}{n-2}\sum_{i=1}^n(y_i-\hat\alpha-\hat\beta x_i)^2$$
On simplification,
\begin{align}
(n-2)s^2&=\sum_{i=1}^n (y_i-\hat\alpha-\hat\beta x_i)y_i
\\&=\sum_{i=1}^n \left\{(y_i-\bar y)-\hat\beta(x_i-\bar x)\right\}y_i
\\&=\sum_{i=1}^n (y_i-\bar y)y_i-\hat\beta \sum_{i=1}^n (x_i-\bar x)y_i
\\&=\sum_{i=1}^n (y_i-\bar y)^2-\hat\beta\sum_{i=1}^n (x_i-\bar x)(y_i-\bar y)
\\&=s_{yy}-\hat\beta s_{xy}
\\&=s_{yy}-\hat\beta^2 s_{xx}
\end{align}
Taking $\alpha'=\alpha+\beta\bar x$, joint pdf of $Y=(Y_1,\ldots,Y_n)$ for $(y_1,\ldots,y_n)\in\mathbb R^n$ is
$$f_{Y}(y_1,\ldots,y_n)=\frac{1}{(\sigma\sqrt{2\pi})^n}\exp\left[-\frac{1}{2\sigma^2}\sum_{i=1}^n\left(y_i-\alpha'-\beta(x_i-\bar x)\right)^2\right]$$
Consider the orthogonal transformation $(y_1,\ldots,y_n)\to(z_1,\ldots,z_n)$ such that
$$\begin{pmatrix}z_1\\z_2\\\vdots\\z_n\end{pmatrix}=\mathbf Q\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}\,,$$
where $$\mathbf Q=\left[\begin{matrix}\frac{1}{\sqrt n}&\frac{1}{\sqrt n}&\cdots&\frac{1}{\sqrt n}\\\frac{x_1-\bar x}{\sqrt{s_{xx}}}&\frac{x_2-\bar x}{\sqrt{s_{xx}}}&\cdots&\frac{x_n-\bar x}{\sqrt{s_{xx}}}\\\vdots&\vdots&\cdots&\vdots\end{matrix}\right]$$ is an $n\times n$ orthogonal matrix with its first two rows fixed.
Then, $$z_1=\frac{1}{\sqrt n}\sum_{i=1}^n y_i=\sqrt{n}\bar y$$ and $$z_2=\frac{\sum (x_i-\bar x)y_i}{\sqrt{s_{xx}}}=\frac{s_{xy}}{\sqrt{s_{xx}}}=\hat\beta\sqrt{s_{xx}}$$
Note that $\sum\limits_{i=1}^n y_i^2=\sum\limits_{i=1}^n z_i^2$ by virtue of orthogonal transformation, which leads to
\begin{align}
\sum_{i=1}^n(y_i-\alpha'-\beta(x_i-\bar x))^2&=\sum_{i=1}^n y_i^2+n\alpha'^2+\beta^2\sum_{i=1}^n(x_i-\bar x)^2-2\alpha'n\bar y-2\beta\sum_{i=1}^n(x_i-\bar x)y_i
\\&=\sum_{i=1}^n z_i^2 +n\alpha'^2+\beta^2 s_{xx}-2\alpha'\sqrt n z_1-2\beta z_2\sqrt{s_{xx}}
\\&=(z_1-\sqrt n\alpha')^2+(z_2-\beta\sqrt{s_{xx}})^2+\sum_{i=3}^nz_i^2
\end{align}
For $(z_1,\ldots,z_n)\in\mathbb R^n$, joint density of $Z=(Z_1,\ldots,Z_n)$ becomes
$$f_{Z}(z_1,\ldots,z_n)=\frac{1}{(\sigma\sqrt{2\pi})^n}\exp\left[-\frac{1}{2\sigma^2}\left\{(z_1-\sqrt n\alpha')^2+(z_2-\beta\sqrt{s_{xx}})^2+\sum_{i=3}^nz_i^2\right\}\right]\,,$$
so that $Z_1,Z_2,\ldots,Z_n$ are independently distributed with
\begin{align}
Z_1&\sim\mathcal N(\sqrt n\alpha',\sigma^2)
\\Z_2&\sim\mathcal N(\beta\sqrt{s_{xx}},\sigma^2)
\\Z_i&\sim\mathcal N(0,\sigma^2)\qquad,\,i=3,4,\ldots,n
\end{align}
Now,
\begin{align}
(n-2)s^2&=s_{yy}-\hat\beta^2s_{xx}
\\&=\sum_{i=1}^ny_i^2-n\bar y^2-\hat\beta^2s_{xx}
\\&=\sum_{i=1}^nz_i^2-z_1^2-z_2^2
\\&=\sum_{i=3}^nz_i^2
\end{align}
And you have that $Z_3,\ldots,Z_n\sim\mathcal N(0,\sigma^2)$, independently.
This implies $$\sum_{i=3}^n\frac{Z_i^2}{\sigma^2}\sim\chi^2_{n-2}$$
Or in other words, $$\frac{(n-2)s^2}{\sigma^2}\sim\chi^2_{n-2}$$
