Nitsche's method is related to discontinuous Galerkin methods (indeed, as Wolfgang points out, it is a precursor to these methods), and can be derived in a similar fashion. Let's consider the simplest problem, Poisson's equation:
$$
\left\{\begin{aligned}
-\Delta u &= f \qquad\text{on }\Omega,\\
u &= g \qquad\text{on }\partial\Omega.
\end{aligned}\right.
\tag{1}
$$
We are now looking for a variational formulation that
- is satisfied by the (weak) solution $u\in H^1(\Omega)$ (i.e., consistent),
- is symmetric in $u$ and $v$,
- admits a unique solution (which means that the bilinear form is coercive).
We start as usual by taking the strong form of the differential equation, multiplying by a test function $v\in H^1(\Omega)$ and integrating by parts. Starting with the right-hand side, we obtain
$$
\begin{aligned}
(f ,v) = (-\Delta u,v)&=(\nabla u,\nabla v) - \int_{\partial\Omega} \partial_\nu u v\,ds \\
&= (\nabla u,\nabla v) - \int_{\partial\Omega} \partial_\nu u v\,ds - \int_{\partial\Omega} (u-g)\partial_\nu v\,ds
\end{aligned}
$$
where in the last equation we have added the productive zero $0=u-g$ on the boundary. Rearranging the terms to separate linear and bilinear forms now gives a variational equation for a symmetric bilinear form that is satisfied for the solution $u\in H^1(\Omega)$ of $(1)$.
The bilinear form is however not coercive, since you cannot bound it from below for $u=v$ by $c\|v\|_{H^1}^2$ (as we don't have any boundary conditions for arbitrary $v\in H^1(\Omega)$, we cannot use Poincaré's inequality as usual -- this means we can make the $L^2$ part of the norm arbitrarily large without changing the bilinear form). So we need to add another (symmetric) term that vanishes for the true solution: $\eta\int_{\partial\Omega} (u-g)v\,ds$ for some $\eta>0$ large enough. This leads to the (symmetric, consistent, coercive) weak formulation: Find $u\in H^1(\Omega)$ such that
$$
(\nabla u,\nabla v) - \int_{\partial\Omega} \partial_\nu u v\,ds - \int_{\partial\Omega} u\partial_\nu v\,ds +\eta\int_{\partial\Omega} u v\,ds = -\int_{\partial\Omega} g\partial_\nu v\,ds + \eta\int_{\partial\Omega} g v\,ds + \int_\Omega f v\,dx
\qquad\text{for all }v\in H^1(\Omega).
$$
Taking instead of $u,v\in H^1(\Omega)$ discrete approximations $u_h,v_h\in V_h\subset H^1(\Omega)$ yields the usual Galerkin approximation. Note that since it's non-conforming due to the boundary conditions (we are looking for the discrete solution in a space that is larger than the one we sought the continuous solution in), one cannot deduce well-posedness of the discrete problem from that of the continuous problem. Nitsche now showed that if $\eta$ is chosen as $ch^{-1}$ for $c>0$ sufficiently large, the discrete problem is in fact stable (with respect to a suitable mesh-dependent norm).
(This is not Nitsche's original derivation, which predates discontinuous Galerkin methods and starts from an equivalent minimization problem. In fact, his original paper does not mention the corresponding bilinear form at all, but you can find it in, e.g., Freund and Stenberg, On weakly imposed boundary conditions for second-order problems, Proceedings of the Ninth Int. Conf. Finite Elements in Fluids, Venice 1995. M. Morandi Cecchi et al., Eds. pp. 327-336.)
