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.
