This is easy to do with [PARI/GP](https://pari.math.u-bordeaux.fr/).
Here is my code
```
p(n) = binomial(n+2,2);
Y(k,b) = sumnum(l=0, (2*l+3)/(p(k)+p(l))^b);
Z(a,b) = sumnum(k=0, (2*k+3)/p(k)^a*Y(k,b));
default(realprecision,57);
print(2*Z(2,2)+4*Z(1,3))
/* 16.0000000000000000000000000000000000000000000000000000000 */
```
It takes four seconds for 57 decimal places. The
[documentation](https://pari.math.u-bordeaux.fr/doc.html)
has some information about the methods being used and there
are several summation functions other than `sumnum`. For
example, I originally used `sumpos` but `sumnum` is better.