This integral is divergent. How to use NIntegrate to see how it grows?

I need to explore the behavior as ϵ goes to 0

I think the behavior is that of $\ln(\pi-x)$ as $x$ approaches $\pi$.

We only need to look at one term.

expr = 1/((3 + Cos[x]) (Sqrt[(3 + Cos[x])^2 - 4]));
Apart[expr]

It is only the first term which makes the integral improper at $x=\pi$. All others are ok.

Integrating the first term, but using series form to make it more clear to see (we only want to see what happens at $\pi$)

Normal@Series[Sqrt[5 + 6 Cos[x] + Cos[x]^2]/(8 (1 + Cos[x])), {x, Pi, 6}]

res=Integrate[%, x]

We see now why it blows up at $x=\pi$ due to the $\ln(\pi-x)$ in the end there. When $x=\pi$ it becomes $\ln(0)$.

 Limit[res, x -> Pi]

So the behavior is simply that of $\ln(\pi-x)$ as $x$ approaches $\pi$?

We can calculate an indefinite integral $\int_a f(x) \; dx$ with NDSolve[{y'[x] == f[x], y[a] == 0}, y, {x, a, b}] yields the integral from a to x, up to x <= b or until NDSolve runs into a singularity or stiffness.

{sol} = NDSolve[{y'[x] == 
    2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])), y[Pi - 0.3] == 0}, 
  y, {x, Pi - 0.3, Pi}];

NDSolve::ndsz: At x == 3.141592536731091`, step size is effectively zero; singularity or stiff system suspected.

Quick visualization of the steps:

ListPlot[y /. sol, PlotRange -> All]

To make a table, let's examine some of the data stored in sol.

dom = y["Domain"] /. sol    (* domain of the solution *)
(*  {{2.8415926535897933`, 3.141592536731091`}}  *)

xmax = dom[[1, 2]];         (* upper limit of the domain *)

Log10[Pi - xmax]            (* just to show the power of 10 *)
(*  -6.93234  *)

This creates a table of the log of the change in x from Pi and the value of the integral:

tab = Table[{N@k, N[y[Pi - 10^-k] /. sol]}, {k, Range[-Ceiling@Log10[Pi - xmax]]}]
(*
  {{1., 0.769312}, {2., 2.39654}, {3., 4.02471},
   {4., 5.65288}, {5., 7.28106}, {6., 8.90923}}
*)

ListPlot[tab, Frame -> True, FrameLabel -> {HoldForm[Log10[Δx]], HoldForm[y[π - Δx]]}]

Repeating the process with a higher working precision extends the domain closer to Pi. (The step size of the integration process is limited by the precision of the numbers used.)

{sol} = NDSolve[{y'[x] == 
    2/((3 + Cos[x])*(Sqrt[((3 + Cos[x])^2 - 4)])), y[Pi - 0.3] == 0}, 
  y, {x, Pi - 0.3, Pi},
  PrecisionGoal -> 16, WorkingPrecision -> 50]

The exact solution can be obtained in such a way:

j = Integrate[ 1/((3 + Cos[x])*Sqrt[(3 + Cos[x])^2 - 4]), {x, Pi/2, Pi - epsilon},
Assumptions -> epsilon > 0 && epsilon < 1/2];
k = Series[j, {epsilon, 0, 2}, Assumptions -> epsilon > 0 && epsilon < 1/2];
FullSimplify[k]

$$O\left(\text{epsilon}^3\right) +\frac{13 \text{epsilon}^2}{192 \sqrt{2}}+\frac{\log \left(\frac{1}{243} (-32) \left(8 \sqrt{10}-37\right)\right)-2 \log (\text{epsilon})}{4 \sqrt{2}} $$