I understand that if a pair quark/antiquark, out of the vacuum, tries to separate then the energy increases, another pair is produced, and we finish with two mesons or generically two hadron jets. But in numerical simulations we can just put a single quark and calculate the field, can we? What does it happen, does the total energy diverge? Or does it increases up to, say, some "constituient" mass?
The total energy diverges, but you can't simulate a single static colored source because this conflicts with periodic boundary conditions. You need to have the box color neutral. The easiest way to do this is to put a color source and antisource far apart, and this was the situation in the earliest lattice QCD simulations.
The color fields in these simulations run in a narrow tube from the source to the antisource. The energy diverges as a linear function of the separation, and the coefficient of the divergence is the QCD string tension. In a simulation with dynamical quarks, quarks will pop out of the vacuum to screen the two sources, so that they will be noninteracting at long distances.
As usual with these things, the presence of QFT and a non-Abelian gauge theory makes life hard, so let's take the prototypical theory that we can actually calculate with easily: Maxwell's equations.
The question is then something like "what happens if I put a single electron in a (classical) lattice simulation?". Immediately, one realises that this question is nonsense because Maxwell's equations only define the time evolution of physical states, which in this case implies things that satisfy Gauss' law (which is a constaint!): $$\nabla \cdot E = \rho.$$ In other words, in a simulation one must start with a physical state before one can run it. The Gauss constraint can be solved in this case, and gives the obvious $1/r^2$ decaying $E$ field.
At this point, we can zoom back out to QCD, and think about what changes. Solving for the Gauss constraint (i.e. finding an element of the physical Hilbert space) is now highly non-trivial, but heuristically we can have a guess at what might go wrong. The basic problem (as Scott Carnahan points out in the comments) is anti-screening. The colour field self-interacts and gets stronger at larger $r$ --- so the total energy diverges, which does not correspond to any element of the physical Hilbert space.
Incidentally, this general problem of needing to solve a constraint equation before being able to start a simulation is quite hard, and numerically one has to be quite clever to avoid instabilities --- many time evolution codes are only stable if you start them with things that satisfy the gauge constraints. Apart from gauge field simulations, simulations of GR also have this problem.