I can understand that if Y1~Yn are random samples from N(μ,σ), then the sum of squares of difference between Yi and bar(Y) divided by sigma^2 follows chi-square distribution with n-1 degress of freedom.
But I can't easily prove that sum of squares of error follows chi-square distribution with n-2 degrees of freedom because it is difference between Yi and estimated Yi.
How can i prove that not using matrix form?
Without involving asymptotic results, you have to assume that the error terms follow Gaussian (uncorrelated) distribution with a constant variance, i.e., $\epsilon \sim N(0, \sigma^2)$. Next, you can show that each fitted value $\hat{Y}_i$ is normally distributed as well, which is follows from $$ \hat{Y}_i = Y_i + \sum_{i=1}^n w_i \epsilon _i, \quad i=1,...n, $$
thus $\hat{Y}_i \sim N(Y_i, \sigma^2 \sum_{i=1}^nw_i^2)$, for each $i$ and they are uncorrelated as well. Now, recall that if $X_i \sim N(0, \sigma^2)$, i.i.d then $$ \sum_{i=1}^n\frac{(X_i - \mu)^2}{\sigma^2} \sim \chi^2(n). $$ Thus, $$ \sum_{i=1}^n \frac{(Y_i - \hat{Y_i})^2}{\sigma^2}\sim \chi^2{(n-2)}. $$ The "reduction" of $2$ df stems from the estimation of $\beta_0$ and $\beta_1$. So, strictly speaking the residuals sum of squares follows $\sigma^2 \chi^2_{(n-2)}$ distribution. To figure out the exact form of $w_i$ you will need some simple algebra of the OLS estimators, but note that $w_i = g(x_i, \bar{x}_n, \sum (x_i -\bar{x}_n)^2)$ so they are considered "constant" terms.