Let lowercase bold letters denote vectors. The linear model is:
$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$
In matrix notation for all observations $i=1, \ldots, n$:
$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$
The OLS estimator for $\mathbf{b}$ is:
$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$
Substituting:
$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$
What's the variance of our estimator $\hat{\mathbf{b}}$? Take the variance (conditional on X)$:
\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}
We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon}$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.
Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}
$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.
You can estimate $\sigma^2$ using the residuals from your OLS regression. Let $\hat{e}_i$ be the residual (as opposed to error-term $\epsilon_i$) for observation $i$. That is:
$$ \hat{e}_i = y_i - \mathbf{x}_i' \hat{\mathbf{b}} $$ Note that this is based upon estimate $\hat{\mathbf{b}}$ instead of true values $\mathbf{b}$. Then the usual estimator for $\sigma^2$ is: $$\hat{\sigma}^2 = \frac{1}{n-k} \sum_i \hat{e}_i^2$$ where $k$ is your number of regressors (including a constant).