Following OP comment requesting more details to Maximilian Matthé's answer:
Just like mentioned, the key to your question is the Bilinear Transform, where you substitute $s$ (or $j\omega$ in your case) with $\frac{2}{T}\frac{z-1}{z+1}$
This transform will have for effect to convert your transfer function from the continuous $s$ domain to the discrete $z$ domain.
I assume that you are fimiliar with the $z$ transform and that you know that a delay in the time domain is equivalent to multiplying by $z^{-1}$ in the z domain. $Z\{f[x-n]\} = Z\{f[x]\}z^{-n}$
The implementation of a discrete filter (FIR or IIR) will happen in the time domain. Your software will implement something like
$$y[n] = b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]...+b_{N}x[n-N]-a_{1}y[n-1] - a_{2}y[n-2]...-a_{N}y[n-N]$$
In my previous equations :
- $x$ is the input
- $y$ is the output
- $n$ is the index of the actual value in your buffer. $n-1$ being the previous value and so on.
If you bring that equation in the $z$ domain, you will have an equivalent like that.
$$Y = b_{0}X + b_{1}Xz^{-1}+b_{2}Xz^{-2}...+b_{N}Xz^{-N} -a_{1}Yz^{-1}-a_{2}Yz^{-2}...-a_{N}Yz^{-N}$$
After some algebra, you'll get the following form, which correspond to a transfer function (ratio of output on input)
$$\frac{Y}{X} = \frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}...+b_{N}z^{-N}}{1+a_{1}z^{-1}+a_{2}z^{-2}...+a_{N}z^{-N}}$$
This transfer function is what Maximilian Matthé answer refers to. In it's condensed form :
$$\frac{Y}{X}=\frac{\sum_{i=0}^{N}b_{i}z^{-i}}{1+\sum_{i=1}^{N}a_{i}z^{-i}}$$
Designing a discrete filter correspond to determining the coefficients of that transfer function. e.g. Finding the values of $a$ and $b$.
The value of your coefficients will change according to your sampling frequency, and you can confirm that considering that the Bilinear Transform introduce the sampling time $T$ in your equation.
You might want to consider the fact that the Bilinear transform is an approximation. If you look at Wikipedia's page, you will find this :

To compensate this approximation, you can substitute $\omega$ with
$$\omega_d = \frac{2}{T}\tan^{-1}(\frac{2\omega}{T})$$
This substitution is called frequency warping and will have a slight effect on the frequency response of your filter, but will most likely be negligible.
Finally, if you have access to matlab and the signal processing toolbox, you might want to have a look to bilinear() which will do all of this in a single command.