Maple does it correctly

int(sqrt((cos(t1) - cos(t2))^2 + (sin(t1) - sin(t2))^2), t1 = 0 .. 2*Pi, t2 = 0 .. 2*Pi)/(2*Pi)^2;
                           4 
                           --
                           Pi

The FullSimplify expression has the same shape but gives the correct result

In[5]:= FullSimplify @ 
Sqrt[(Cos[\[Theta]1] - Cos[\[Theta]2])^2 + (Sin[\[Theta]1] - 
Sin[\[Theta]2])^2]

Out[5]= Sqrt[2 - 2 Cos[\[Theta]1 - \[Theta]2]]

In[6]:= Plot3D[Sqrt[
2 - 2 Cos[\[Theta]1 - \[Theta]2]], {\[Theta]1, 0, 
2 \[Pi]}, {\[Theta]2, 0, 2 \[Pi]}]

You are right, the result is wrong. Here are three ways of calculating the integral, two give the right answer and one does not

f = Sqrt[(Cos[\[Theta]1] - Cos[\[Theta]2])^2 + (Sin[\[Theta]1] - 
       Sin[\[Theta]2])^2];
Integrate[
 FullSimplify[f], {\[Theta]1, 0, 2 \[Pi]}, {\[Theta]2, 0, 2 \[Pi]}]
Integrate[f, {\[Theta]1, \[Theta]2} \[Element] 
  Rectangle[{0, 0}, {2 Pi, 2 Pi}]]
Integrate[f, {\[Theta]1, 0, 2 \[Pi]}, {\[Theta]2, 0, 2 \[Pi]}]

However, it is a sad fact that symbolic integration is prone to erros due to branch cuts. The function we are integrating is non-analytical on the diagonal, and I suppose this is the root of the problem. The primitive with respect to one variable is discontinuous on the diagonal:

Integrate[f, \[Theta]1]
Plot3D[%, {\[Theta]1, 0, 2 \[Pi]}, {\[Theta]2, 0, 2 \[Pi]}]