Different regions of your integrand present different challenges to the integrator.

As a first simple experiment I subdivided the integrand into the interior, the six "faces" and the eight "corners." Some of those regions appear to integrate relatively nicely while others seem to have substantially more trouble. Actually I'm not certain that there might need to be further subdivisions near t=Pi/2.

Here is what I've been looking at:

In[1]:= d=1/30;

f=Simplify[2*1/(3!*π^2)*r*(2 r*Sin[t]*(r^2+z1^2-2*r*z1*Cos[t]))/((1+r^2 z1^2-2 r z1 Cos[t])*

Sqrt[((1+r^2*Sin[t]^2-r^2*Cos[t]^2)^2+4*r^4*Sin[t]^2*Cos[t]^2)]*Sqrt[1-z1^2]*(1-r^2))]

Interior

In[3]:= NIntegrate[f,{z1,-1+d,1-d},{r,0+d,1-d},{t,d Pi,π-d Pi}]

Faces

In[4]:= NIntegrate[f,{z1,-1,-1+d},{r,d,1-d},{t,d Pi,π-d Pi}]

In[5]:= NIntegrate[f,{z1,1-d,1},{r,d,1-d},{t,d Pi,π-d Pi}]

In[6]:= NIntegrate[f,{z1,-1+d,1-d},{r,0,d},{t,d Pi,π-d Pi}]

In[7]:= NIntegrate[f,{z1,1-d,1},{r,1-d,1},{t,d Pi,π-d Pi}]

In[8]:= NIntegrate[f,{z1,-1+d,1-d},{r,d,1-d},{t,0,d Pi}]

In[9]:= NIntegrate[f,{z1,1-d,1},{r,d,1-d},{t,π-d Pi,Pi}]

Corners

In[10]:= NIntegrate[f,{z1,-1,-1+d},{r,0,d},{t,0,d Pi}]

In[11]:= NIntegrate[f,{z1,-1,-1+d},{r,0,d},{t,Pi-d Pi,Pi}]

In[12]:= NIntegrate[f,{z1,-1,-1+d},{r,1-d,1},{t,0,d Pi}]

In[13]:= NIntegrate[f,{z1,-1,-1+d},{r,1-d,1},{t,Pi-d Pi,Pi}]

In[14]:= NIntegrate[f,{z1,1-d,1},{r,0,d},{t,0,d Pi}]

In[15]:= NIntegrate[f,{z1,1-d,1},{r,0,d},{t,Pi-d Pi,Pi}]

In[16]:= NIntegrate[f,{z1,1-d,1},{r,1-d,1},{t,0,d Pi}]

In[17]:= NIntegrate[f,{z1,1-d,1},{r,1-d,1},{t,Pi-d Pi,Pi}]

Inspect the denominator of the simplified expression carefully and look for regions where the denominator will go to zero. I believe those will tend to pose the greatest challenges to the integrator.

As you increase the PrecisionGoal this will further expose which of these regions have increased problems.