If d Cos[u] is not accepted as integration-variable you should change it to ( - Sin[u] du ). You have to change the boundaries as well
Then: your inner integrand is really complicated containing transcedent functions.
So I tried it numerically in total:
f1[R_, r_?NumericQ] :=
NIntegrate[(R - r*Cos[u])*(3 Cos[u]^2 - 1)*
E^(-1.5566*(r^2 + R^2 - 2*r*R*Cos[u])^0.5) (-Sin[u]), {u, -Pi, 0}]
and
NIntegrate[(Sqrt[15]/4)*f1[3.9, r]*r^4*E^(-3.313066*r), {r, 0, 3.9}]
