My impression is that the prime numbers in $S(3k+1)$ are exactly in the proportion they should be. Note that numbers in $S(3k+1)$ are smaller than $10^{D(3k+1)-1}$ because they have $D(3k+1)$ digits.

By the prime number theorem, primes up to $10^{D(3k+1)-1}$ are in proportion $1/(D(3k+1)-1)\log(10) $. In the subset $S(3k+1)$, which seems "arithmetically random", I would expect the same proportion, that is $$|P(3k+1)| \sim \frac{|S(3k+1)|}{(D(3k+1)-1)\log(10)} \ \ \ \ :NO:$$ But hey (thanks Noam D. Elkies), we have to remember such numbers are never multiples of $3$. This is the only "easy" constraint we have. Taking into account this fact, there is a factor $3/2$ to take into account $$ |P(3k+1)| \sim \frac{3|S(3k+1)|}{2(D(3k+1)-1)\log(10)}$$ The factor is explained in an edit at the end. Using the fact that $S(3k+1) \sim (3k+1)!$, we get the conjecture $$|P(3k+1)| \sim \frac{3(3k+1)!}{2(D(3k+1)-1)\log(10)}$$ I think this is the "naive guess", but actually it is not consistent with yours.

Indeed... If I understand the notation correctly, there is an even worse problem in your limit! Let's estimate $D(n)$. Since a number $\ell \le n$ has about $\log_{10}(\ell)$ digits, we have $$D(n) \sim \sum_{\ell = 0}^{n}\log_{10}(\ell) = \frac{ \log ( n!)}{ \log(10) } $$ Using the Stirling approximation $$ D(n) \sim \frac{ \log(n!) }{ \log(10) } \sim \frac{n \log(n) }{ \log(10)} = n \log_{10}(n)$$

This means for $\log_{10}(n)$ big enough, $D(n) \ge 2n$. You don't see this phenomenon for small $n$ since $\log_{10}(n)$ is at most 1 in your case. But hey... this means $$(D(3k+1) -1)! \ge (6k+1)! >> (3k+1)! = |S(3k+1)| $$ In other words, the estimate you are trying to give (if I understood notation correctly) is even bigger than the set you extract $P(3k+1)$ from.

Final remark. Notice that we can plug the Stirling approximation for $D(n)$ into the estimate above to get $$ P(3k+1) \sim \frac{ 3(3k+1)!}{ 2(3k+1) \log(3k+1) \log(10) } = \frac{ 3(3k)! }{ 2 \log(3k+1) \log(10) } $$

However, this performs much worse than your estimate for the first three values, getting $3,240,102664$. But I would not discard an asymptotic estimate because it doesn't work on the first values! It is possible though that I have been too sloppy :)

EDIT The factor $3/2$ can be explained in this way. Call $A(j)_{\le x}$ the subset of numbers congruent to $j$ modulo $3$ smaller than $x$, for $j=0,1,2$. If we have $m$ primes in $A(0)_{\le 3x} \sqcup A(1)_{\le 3x} \sqcup A(2)_{\le 3x} $, then the primes up to $3x$ are $m/3x$, while the primes up to $3x$ in $A(1)_{\le 3n} \sqcup A(2)_{\le 3n}$ are $$\frac{m}{2x} = \frac{3}{2} \cdot \frac{m}{3x}$$ Thus the density in $A(1)_{\le 3x} \sqcup A(2)_{\le x}$ is (as the intuition suggests) higher, and there is precisely a factor of $3/2$. Also, by the extended prime number theorem, $A(1)_{\le 3x}$ and $A(2)_{\le 3x}$ have approximately the same number of primes. This means $A(1)_{\le 3x}$ contains about $m/2$ primes, resulting in a density of $$\frac{m/2}{x} = \frac{m}{2x} = \frac{3}{2} \frac{m}{3x}$$ Explaining the factor.

My impression is that the prime numbers in $S(3k+1)$ are exactly in the proportion they should be. Note that numbers in $S(3k+1)$ are smaller than $10^{D(3k+1)-1}$ because they have $D(3k+1)$ digits.

By the prime number theorem, primes up to $10^{D(3k+1)-1}$ are in proportion $1/(D(3k+1)-1)\log(10) $. In the subset $S(3k+1)$, which seems "arithmetically random", I would expect the same proportion, that is $$|P(3k+1)| \sim \frac{|S(3k+1)|}{(D(3k+1)-1)\log(10)} \ \ \ \ :NO:$$ But hey (thanks Noam D. Elkies), we have to remember such numbers are never multiples of $3$. This is the only "easy" constraint we have. Taking into account this fact, there is a factor $3/2$ to take into account $$ |P(3k+1)| \sim \frac{3|S(3k+1)|}{2(D(3k+1)-1)\log(10)}$$ The factor is explained in an edit at the end.

(EDITED) We are left with estimating $S(3k+1)$. We have $D(3k+1)$ digits, and it's not hard to see that are almost equally distributed among the ten possible digits. There is indeed a bias, but I will neglect it (see below). The permutations of $D(3k+1)$ things divided in 10 equal groups is $$\frac{ D(3k+1)!}{ \left ( \frac{ D(3k+1)}{10} \right ) !^{10} } $$ Using the complete stirling approximation, after a few calculations we get $$ \frac{10^{D(3k+1)+5}}{ (2 \pi D)^{9/2} )}$$ Since about $1/10$ of the numbers will start with $0$, we get in the end $$ S(3k+1) \sim \frac{10^{D(3k+1)+4}}{ (2 \pi D)^{9/2} }$$ Let's get back to estimate $P(3k+1)$. We have $$|P(3k+1)| \sim \frac{3}{2} \cdot \frac{10^{D(3k+1)+4}}{ (2 \pi D(3k+1) )^{9/2} ( D(3k+1) -1) \log(10)}$$ I think this is the "naive guess", but actually it is not consistent with yours.

Now let's estimate $D(n)$. Since a number $\ell \le n$ has about $\log_{10}(\ell)$ digits, we have $$D(n) \sim \sum_{\ell = 0}^{n}\log_{10}(\ell) = \frac{ \log ( n!)}{ \log(10) } $$ Using the Stirling approximation $$ D(n) \sim \frac{ \log(n!) }{ \log(10) } \sim \frac{n \log(n) }{ \log(10)} = n \log_{10}(n)$$ If we plug the Stirling approximation for $D(n)$ into the estimate above we get $$|P(m)| \sim \frac{10^{m \log_{10}(m) +4}}{ (2 \pi)^{9/2} (m \log_{10}(m) )^{11/2} \log(10)}$$ $$|P(m)| \sim \frac{ 3 \cdot 10000}{ 2 \cdot ( 2 \pi)^{9/2} \cdot \log(10) }\frac{m^{m-\frac{11}{2} }}{ \log_{10}(m)^{11/2}}$$ $$|P(m)| \sim 1.667 \cdot \frac{m^{m-\frac{11}{2} }}{ \log_{10}(m)^{11/2}}$$

I would be curious to see how this performs!

EDIT The factor $3/2$ can be explained in this way. Call $A(j)_{\le x}$ the subset of numbers congruent to $j$ modulo $3$ smaller than $x$, for $j=0,1,2$. If we have $m$ primes in $A(0)_{\le 3x} \sqcup A(1)_{\le 3x} \sqcup A(2)_{\le 3x} $, then the primes up to $3x$ are $m/3x$, while the primes up to $3x$ in $A(1)_{\le 3n} \sqcup A(2)_{\le 3n}$ are $$\frac{m}{2x} = \frac{3}{2} \cdot \frac{m}{3x}$$ Thus the density in $A(1)_{\le 3x} \sqcup A(2)_{\le x}$ is (as the intuition suggests) higher, and there is precisely a factor of $3/2$. Also, by the extended prime number theorem, $A(1)_{\le 3x}$ and $A(2)_{\le 3x}$ have approximately the same number of primes. This means $A(1)_{\le 3x}$ contains about $m/2$ primes, resulting in a density of $$\frac{m/2}{x} = \frac{m}{2x} = \frac{3}{2} \frac{m}{3x}$$ Explaining the factor.

BIAS. If $(3k+1)$ is for example $2002$, one thousand the numbers will start with $1$, so there will be a relative abundance of $1$. Here $D(2002) \sim 8000$ so each digit should appear $800$ times, but instead the digit $1$ appears $1000$ of times + $D(999)/10 \simeq 400$, while the other digits appear a little less than $800$ times.

My impression is that the prime numbers in $S(3k+1)$ are exactly in the proportion they should be. Note that numbers in $S(3k+1)$ are smaller than $10^{D(3k+1)-1}$ because they have $D(3k+1)$ digits.

By the prime number theorem, primes up to $10^{D(3k+1)-1}$ are in proportion $1/(D(3k+1)-1)\log(10) $. In the subset $S(3k+1)$, which seems "arithmetically random", I would expect the same proportion, that is $$|P(3k+1)| \sim \frac{|S(3k+1)|}{(D(3k+1)-1)\log(10)}$$ Using the fact that $S(3k+1) \sim (3k+1)!$, we get the conjecture $$|P(3k+1)| \sim \frac{(3k+1)!}{(D(3k+1)-1)\log(10)}$$ I think this is the "naive guess", but actually it is not consistent with yours.

Indeed... If I understand the notation correctly, there is an even worse problem in your limit! Let's estimate $D(n)$. Since a number $\ell \le n$ has about $\log_{10}(\ell)$ digits, we have $$D(n) \sim \sum_{\ell = 0}^{n}\log_{10}(\ell) = \frac{ \log ( n!)}{ \log(10) } $$ Using the Stirling approximation $$ D(n) \sim \frac{ \log(n!) }{ \log(10) } \sim \frac{n \log(n) }{ \log(10)} = n \log_{10}(n)$$

This means for $\log_{10}(n)$ big enough, $D(n) \ge 2n$. You don't see this phenomenon for small $n$ since $\log_{10}(n)$ is at most 1 in your case. But hey... this means $$(D(3k+1) -1)! \ge (6k+1)! >> (3k+1)! = |S(3k+1)| $$ In other words, the estimate you are trying to give (if I understood notation correctly) is even bigger than the set you extract $P(3k+1)$ from.

Final remark. Notice that we can plug the Stirling approximation for $D(n)$ into the estimate above to get $$ P(3k+1) \sim \frac{ (3k+1)!}{ (3k+1) \log(3k+1) \log(10) } = \frac{ (3k)! }{ \log(3k+1) \log(10) } $$

However, this performs much worse than your estimate for the first three values, getting $2,160,68443$. But I would not discard an asymptotic estimate because it doesn't work on the first values! It is possible though that I have been too sloppy :)

My impression is that the prime numbers in $S(3k+1)$ are exactly in the proportion they should be. Note that numbers in $S(3k+1)$ are smaller than $10^{D(3k+1)-1}$ because they have $D(3k+1)$ digits.

By the prime number theorem, primes up to $10^{D(3k+1)-1}$ are in proportion $1/(D(3k+1)-1)\log(10) $. In the subset $S(3k+1)$, which seems "arithmetically random", I would expect the same proportion, that is $$|P(3k+1)| \sim \frac{|S(3k+1)|}{(D(3k+1)-1)\log(10)} \ \ \ \ :NO:$$ But hey (thanks Noam D. Elkies), we have to remember such numbers are never multiples of $3$. This is the only "easy" constraint we have. Taking into account this fact, there is a factor $3/2$ to take into account $$ |P(3k+1)| \sim \frac{3|S(3k+1)|}{2(D(3k+1)-1)\log(10)}$$ The factor is explained in an edit at the end. Using the fact that $S(3k+1) \sim (3k+1)!$, we get the conjecture $$|P(3k+1)| \sim \frac{3(3k+1)!}{2(D(3k+1)-1)\log(10)}$$ I think this is the "naive guess", but actually it is not consistent with yours.

Indeed... If I understand the notation correctly, there is an even worse problem in your limit! Let's estimate $D(n)$. Since a number $\ell \le n$ has about $\log_{10}(\ell)$ digits, we have $$D(n) \sim \sum_{\ell = 0}^{n}\log_{10}(\ell) = \frac{ \log ( n!)}{ \log(10) } $$ Using the Stirling approximation $$ D(n) \sim \frac{ \log(n!) }{ \log(10) } \sim \frac{n \log(n) }{ \log(10)} = n \log_{10}(n)$$

This means for $\log_{10}(n)$ big enough, $D(n) \ge 2n$. You don't see this phenomenon for small $n$ since $\log_{10}(n)$ is at most 1 in your case. But hey... this means $$(D(3k+1) -1)! \ge (6k+1)! >> (3k+1)! = |S(3k+1)| $$ In other words, the estimate you are trying to give (if I understood notation correctly) is even bigger than the set you extract $P(3k+1)$ from.

Final remark. Notice that we can plug the Stirling approximation for $D(n)$ into the estimate above to get $$ P(3k+1) \sim \frac{ 3(3k+1)!}{ 2(3k+1) \log(3k+1) \log(10) } = \frac{ 3(3k)! }{ 2 \log(3k+1) \log(10) } $$

However, this performs much worse than your estimate for the first three values, getting $3,240,102664$. But I would not discard an asymptotic estimate because it doesn't work on the first values! It is possible though that I have been too sloppy :)

EDIT The factor $3/2$ can be explained in this way. Call $A(j)_{\le x}$ the subset of numbers congruent to $j$ modulo $3$ smaller than $x$, for $j=0,1,2$. If we have $m$ primes in $A(0)_{\le 3x} \sqcup A(1)_{\le 3x} \sqcup A(2)_{\le 3x} $, then the primes up to $3x$ are $m/3x$, while the primes up to $3x$ in $A(1)_{\le 3n} \sqcup A(2)_{\le 3n}$ are $$\frac{m}{2x} = \frac{3}{2} \cdot \frac{m}{3x}$$ Thus the density in $A(1)_{\le 3x} \sqcup A(2)_{\le x}$ is (as the intuition suggests) higher, and there is precisely a factor of $3/2$. Also, by the extended prime number theorem, $A(1)_{\le 3x}$ and $A(2)_{\le 3x}$ have approximately the same number of primes. This means $A(1)_{\le 3x}$ contains about $m/2$ primes, resulting in a density of $$\frac{m/2}{x} = \frac{m}{2x} = \frac{3}{2} \frac{m}{3x}$$ Explaining the factor.