Computing the twin prime constant with Mathematica

The problem with your code is that N[Product[...]] calls NProduct[...] which turns integers into reals before Prime can evaluate, causing the error: See Details and Options section of documentation of NProduct.

For your calculation, you actually do not need these. It looks like the result converges around $0.660162$:

Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, 10^k}]], {k, 5}]

{0.665138, 0.66033, 0.66017, 0.660162, 0.660162}

Since the product is monotonic and you are using N at the end, it is not necessary to carry out whole multiplication all the way upto infinity: Depending on the interested accuracy, you can go to higher orders as well. For example,

Table[N[Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, 10^k}], 8], {k, 5}]

yields

{0.66513840, 0.66033029, 0.66017020, 0.66016232, 0.66016185}

which reveals that using only first $10^5$ primes, which is all primes upto $1299709$, is sufficient to get a result with error of the order $10^{-6}$.

For the record: Mathematica fails to do the calculation analytically:

Product[(Prime[i] (Prime[i] - 2))/(Prime[i] - 1)^2, {i, 2, Infinity}]

$\prod _{i=2}^{\infty } \frac{\left(p_i-2\right) p_i}{\left(p_i-1\right){}^2}$

As not all primes are known, Mathematica will have a problem to compute the constant exactly. But the sequence appears to converge quickly to some number $0.660162...$. (Soner posted that already.)

So let me give a way to compute the sequence of partial products more efficientlt.

n = 10000;
seq = FoldList[#1 #2 (#2 - 2)/(#2 - 1)^2 &, 1, Prime[Range[2, n]]];
N[seq[[-1]]]
ListLinePlot[{N[seq], ConstantArray[N[seq[[-1]]], Length[seq]]},
 PlotRange -> {1/2, 1},
 PlotStyle -> {Thick, Dashed}
 ]

0.660162