Mathematica makes this fairly simple.  We can define a region 

    region = RegionDifference[Disk[{a, 0}, a], Disk[{0, 0}, b]]
    (* BooleanRegion[#1 && ! #2 &, {Disk[{a, 0}, a], 
      Disk[{0, 0}, b]}] *)
    
You can visualise this quite simply, for specific values of `a` and `b` with the following line

    Block[{a = 3, b = 1}, Region[region]]



and an integrand

    integrand = x^2 + y^2 - 2 a x + 2 b^2 a x/(x^2 + y^2) - b^2;

and simply apply `Integrate`.  I hope you find the result information!
    
    Assuming[0 < b < a, 
     Integrate[integrand, {x, y} ∈ region]]
    (* 1/24 (-6 a^2 b Sqrt[4 a^2 - b^2] - 
       21 b^3 Sqrt[4 a^2 - b^2] - 12 a^4 π + 8 b^4 π - 
       48 a^2 b^2 ArcCos[Sqrt[b/a]/Sqrt[2]] + 
       48 b^4 ArcCot[b/Sqrt[4 a^2 - b^2]] - 20 a^4 ArcCsc[(2 a)/b] - 
       4 b^4 ArcSec[(2 a)/b] + 44 a^4 ArcSin[b/(2 a)] + 
       48 a^2 b^2 ArcSin[b/(2 a)] - 
       48 a^2 b^2 ArcSin[Sqrt[b/a]/Sqrt[2]] + 
       96 a^2 b^2 ArcTan[Sqrt[-1 + (4 a^2)/b^2]] - 
       48 b^4 ArcTan[Sqrt[-1 + (4 a^2)/b^2]] - 
       16 b^4 ArcTan[b/Sqrt[4 a^2 - b^2]]) *)