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]]) *)

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 informative!

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]]) *)