I have attached an example...it took me some time to put it together, and it requires a bit of discussion.

The reason for using Table was to understand the source of the variations shown in the plot at the bottom of the attachment. The big step function should not be there. Also, the small ripples on the leading flat portion of the plot replicate if I "zoom" in on them by limiting the plot range. They should also not be in there. In fact, I had expected the plot to be a smooth continuation of the curved portion of the plot extrapolating back toward d->0...it never gets to zero but to a small value.

The Tables are for each integral separately and then for the difference between the two integrals and then the ratio of the two integrals. A plot of the Table points show smooth curves with an occasional discontinuity that should also not be in there. My hand calculations using the Table data show that the two integrals differ in the 7th or 8th decimal place of the already small values. In real life, there are common coefficients to the integrals that make them much larger, but those are constants that cancel in the ratio and are common to the other comparisons as simple multipliers.

I suspect my technique is faulty but I cannot find any way to improve the results. I am either seeing noise or the results of an error in performing the integrations and comparisons.

Does the attached file make sense?

Attachments:

Ack, it's the AccuracyGoal. Bring it up to 12, say (even 8 seems to be enough for the range in question). Keeping in mind that the discrepancies are around the 6th digit, it's safe to expect that AccuracyGoal=>6 may not be adequate.

Here are some remarks.

(1) I did some further checking. PrecisionGoal also needs tweaking. This is close to smooth, but slow.

fme[d_?NumericQ] := 
 NIntegrate[
   1000 (r1*r2*(d + z2 - z1)*
      d^2/(r1^2 + r2^2 + z1^2 + z2^2 - 2*z1*z2 - 
          2*r1*r2*Cos[t2 - t1] + d*(d + 2*z2 - 2*z1))^1.5), {t1, 0, 
    2 Pi}, {t2, 0, 2 Pi}, {z1, -1, 1}, {z2, -.27, .27}, {r1, 0, 
    Sqrt[1 - z1^2]}, {r2, 0, Sqrt[.27^2 - z2^2]}, PrecisionGoal -> 12,
    AccuracyGoal -> 15]/
  NIntegrate[
   r1*r2*1000, {t1, 0, 2 Pi}, {t2, 0, 2 Pi}, {z1, -1, 
    1}, {z2, -.27, .27}, {r1, 0, Sqrt[1 - z1^2]}, {r2, 0, 
    Sqrt[.27^2 - z2^2]}, PrecisionGoal -> 12, AccuracyGoal -> 15]

Plot[fme[d], {d, 100, 130}, AxesLabel -> {d, F (Normalized)}, 
 PlotLegends -> {"Sphere"}, PlotStyle -> {Thick}, PlotRange -> Full]

(2) Might or might not help to simplify the denominator and also not recompute it. Can even compute it exactly.

c2 = 27/100;
d2 = 16*Pi^2*
  Integrate[
   r1*r2, {z1, 0, 1}, {z2, 0, c2}, {r1, 0, Sqrt[1 - z1^2]}, {r2, 0, 
    Sqrt[c2^2 - z2^2]}]

Also you can use c2 in the numerator computations.

(3) At lower PrecisionGoal, maybe in the range 8-10, it is reasonably fast and not as jumpy but still shows jumps and smooth sections (at least 6 places out though, since the y axis points are all labeled as 1). I suspect numerical issues wherein error estimates cross a threshold and the routines either do further refinement, or do not, depending on which side of the threshold they are on.

(4) As for your other note, if you saw imaginary parts and jumps due to branch cuts in exact definite integrals involving elliptics, I am not surprised. Should you see them again, please politely remind them that they remain the bane of my existence.