I record this here in the manner of usefulness.

The preprint draft http://www.math.u-bordeaux1.fr/~cohen/hardylw.dvi of Cohen gives a method for computing such constants.

The idea is to turn an Euler product into an Euler sum by taking logarithms, transform by Möbius inversion to turn sums over prime powers into logs of $L$-functions at positive integers, and then evaluate the $L$-functions by "standard" methods. To speed the convergence, one can consider the primes up to $A$ (say 30) separately.

In Sections 1 and 3, Cohen considers how to evaluate said $L$-functions (such methods are standard following Hecke/Lavrik, and are implemented by Dokchitser in GP/PARI or Magma), and in sections 4 and 5 he describes the cases of quadratic and cubic fields more fully. Section 2 gives an overview about the manipulations with sums over prime powers.

As a model case to illustrate, consider (Cohen's Section 4) the case of a quadratic Dirichlet character $\chi$. The base objects are the sums $S_m(\chi)=\sum_p \chi(p)/p^m$, and by Möbius inversion these are equal to

$$S_m(\chi)=\sum_{k\ge 1}^\infty {\mu(k)\over k}\log L(km,\chi^k).$$

However, by splitting the $p$-sum at $A$, one in fact by a similar inversion has

$$S_m(\chi)=\sum_{p\lt A}{\chi(p)\over p^m}+\sum_{k\ge 1}^\infty {\mu(k)\over k}\log L_{p\ge A}(km,\chi^k),$$

and this $\log L_{p\ge A}(km,\chi^k)$ is bounded by $O(1/A^{km})$ (Cohen does this in detail in Section 2.1).

To use the above formula, one needs a method to compute $L(km,\chi^k)$, and as was said, for any $L$-function this is standard (by an approximate functional equation with weights given by an incomplete Mellin transform of the $\Gamma$-factors), and takes time proportional to approximately the square root of the conductor of $\chi^k$ (with mild constants dependent on the desired precision, and also the size of $km$).

In order to obtain the Hardy-Littlewood or Bateman-Horn constants from the base objects $S_m(\chi)$, one needs to piece them together in the proper manner, which is largely bookkeeping (see 4.1 and 4.2 of Cohen), involving Dirichlet convolution of arithmetic sequences.

For higher degree polynomials, instead of Dirichlet $L$-functions one will also need Artin $L$-functions, which Cohen points out could be computed by $\zeta$-quotients, though this would be rather inefficient. Instead, a Magma package of Dokchitser allows one to compute with Artin representations rather easily. The time is still proportional to approximately square root of the conductor, though with higher degree $L$-functions the constant factor(s) are not quite so nice.

I am not sure this really answers the question. In this method, one does not "multiply out" the Euler product, except for the small primes (though one could take $A=2$, if treating small primes differently is not pleasing). The bulk contribution is determined by an $L$-function computation, which I suppose could be seen as an "integral" (in terms of a Mellin transform?) if desired.