Dispersion relation near a Dirac point I have very simple energy spectrum,
$$E_{\pm}(k_x,k_y)=\sqrt{2+2\cos(k_x)\cos(k_y)}.$$
The band gap closes at (for example) $k_x=0, k_y=\pi$. Both $E_+$ and $E_-$ are zero there. This is seen in the plot:

Apparently the function is supposed to be linear in momenta, so I expect the second term in the Taylor expansion to be important. However, I can't seem to compute the Taylor expansion because the terms either vanish or become undefined. A Taylor expansion about the point $(a,b)$ of the function $f(k_x,k_y)$ is,
$$f(k_x,k_y)\approx f(a,b)+\partial_{k_x}f(a,b)(k_x-a)+\partial_{k_y}f(a,b)(k_y-b).$$
Okay. They function $E_+$, is clearly zero at that point. So the first term in the Taylor expansion vanishes. The second term is,
$$\partial_{k_x}E_+(0,\pi)k_x=\frac{-\sin(k_x=0)\cos(k_y=\pi)}{\sqrt{2+2\cos(k_x=0)\cos(k_y=\pi)}}k_x.$$
But this is undefined in the point $k_x=0, k_y=\pi$ because the denominator is zero. So what am I doing wrong? 
Is there a nice way to look at the dispersion relation near a Dirac point?
Perhaps my question is more general, how can I look at a function in a point where the function itself is zero?
 A: 
Is there a nice way to look at the dispersion relation near a Dirac point?

It turns out there is. Namely, since Dirac points are equivalent to having a linear dispersion relation, $ E(k) \propto k $.
Warm-up
Consider the simple absolute-value function,
$$
|x| = \sqrt{x^2} =
\begin{cases}
-x & , ~ x < 0 \\
x & , ~ x \geq 0
\end{cases} ~,
$$
and its derivative,
$$
\frac{d|x|}{dx} = \frac{x}{|x|} =
\begin{cases}
-1 & , ~ x < 0 \\
+1 & , ~ x > 0
\end{cases} ~.
$$
Notice that the derivative is not defined at $ x = 0 $. You cannot make a Taylor expansion around zero by naïvely using the derivative. Nonetheless, we can “approximate” the two branches of the absolute-value function around $ x = 0 $ as
$$
|x| "\approx"
\begin{cases}
-x & , ~ x < 0 \\
+x & , ~ x > 0
\end{cases}~.
$$
The nice thing is that Dirac points have an analogous behaviour.
1D dispersion with Dirac point
To get a better idea, let's first look at a simpler dispersion related to your original one:
$$
E(k) = \sqrt{1 + \cos(k)} ~.
$$
The group velocity,
$$
v_k \equiv \frac{\partial E(k)}{\partial k} = -\frac{1}{2} \frac{\sin(k)}{\sqrt{1 + \cos(k)}} ~,
$$
is obviously undefined at $ k = \pi $, as in your original problem; that means, no Taylor expansion around $ k = \pi $ is possible. Yet, there is a way, like the case of the absolute-value function to approximate the branches of the function around $ k = \pi $, as follows:
$$
\begin{align}
\sqrt{1 + \cos(k)}\Big\vert_{k = \pi} &\equiv \sqrt{1 + \cos(\pi + x)}\Big\vert_{x = 0} \\
&= \sqrt{1 + \cos(\pi) -\frac{1}{2} \cos(\pi) \, x^2 + \mathcal{O}(x^4)} \tag{1} \\
&\approx \sqrt{\frac{1}{2} x^2} = \frac{|x|}{\sqrt{2}} =
\frac{1}{\sqrt{2}} \, |k - \pi| ~, \tag{2}
\end{align}
$$
where at (1), we have used the Taylor expansion of $\cos$ around $ \pi $, and in the approximating step (2), we have neglected nonlinear terms.
Therefore, as $ x \equiv k - \pi $, it is easily seen that the group velocity, or the derivative of $ E(k) $ around $ k = \pi $, reads
$$
v \Big\vert_{k = \pi} \approx \frac{1}{\sqrt{2}}
\begin{cases}
-1 & , ~ k < \pi \\
+1 & , ~ k > \pi ~,
\end{cases}
$$
while the dispersion (2) is obviously linear in $k$; a proper Dirac point.
2D dispersion with Dirac point
Armed with the analysis above, we can tackle your original 2D dispersion,
$$
E(k_x, k_y) = \sqrt{2 + 2 \cos(k_x) \cos(k_y)} ~.
$$
At $(k_x = 0, k_y = \pi)$, we can expand $ E(k_x, k_y) $ as 
$$
\begin{align}
\sqrt{2 + 2 \cos(k_x) \cos(k_y)}\Big\vert_{k_x = 0 , k_y = \pi} &\equiv \sqrt{2 + 2 \cos(x) \cos(\pi + y)}\Big\vert_{x = 0, y = 0} \\
&\approx \sqrt{x^2 + y^2} = |k_x \, (k_y - \pi) | ~, \tag{3}
\end{align}
$$
where we have used
$$
\begin{align}
\cos(x) \Big\vert_{x = 0} = 1 -\frac{1}{2} x^2 + \mathcal{O}(x^4) ~, \\
\cos(\pi + y) \Big\vert_{y = 0} = 1 +\frac{1}{2} y^2 + \mathcal{O}(y^4) ~, \\
\cos(x) \cdot \cos(\pi + y) \Big\vert_{x = 0, y = 0} = -1 + \frac{1}{2} (x^2 + y^2) + \mathcal{O}(x^2 y^2) ~,
\end{align} 
$$
and the definitions, $ x \equiv k_x $ and $ y \equiv k_y - \pi $, while we have kept the leading terms in the expansion.
Therefore,
$$
\begin{align}
v_x \Big\vert_{k_x = 0, k_y = \pi} &\equiv \frac{\partial E(\mathbf{k})}{\partial k_x} \Big\vert_{k_x = 0, k_y = \pi}
= \mathrm{sign}(k_x) \cdot | k_y - \pi | \\
v_y \Big\vert_{k_x = 0, k_y = \pi} &\equiv \frac{\partial E(\mathbf{k})}{\partial k_y} \Big\vert_{k_x = 0, k_y = \pi}
= | k_x | \cdot \mathrm{sign}(k_y - \pi) \tag{4}
\end{align} 
$$
In this way, you recover the bi-linear dispersion relations (3), and group velocities (4) which are characteristics of a Dirac point. This scheme is general in that you can apply it to any Dirac point since the dispersion at Dirac points behaves like the absolute-value function. In a mathematical sense, $ E(k) $ is not smooth at certain points of its domain. Near such nonsmooth points, care must be taken with derivatives, since they are singular or undefined. In such a case, one should “extract the nonsmooth part” or singularity, for instance, by the method I presented above.
As a side remark, notice that a point where the group velocity, $ v_{\mathbf{k}} \equiv \nabla_\mathbf{k} E(\mathbf{k}) $, is undefined, like $(k_x = 0, k_y = \pi)$ here, is called a van-Hove singularity in band-structure theory.
