A likert scale, as the term is typically used, is just an ordinal rating scale. The phrase is often used for a single rating, which might have been called a likert item. Traditionally, the idea was that you would have a set of likert items that all measure the same thing and have the same measurement properties. The result is that you could sum (or average) the items and end up with a good measure of something that approximated a continuous, interval scale (see: levels of measurement).
On the other hand, to simulate data, you need to know the distribution that the data should have. More generally, for simulation studies people generally want to have a data generating process for the resulting distribution. A likert scale is a type of data gathering instrument, not a distribution and not a data generating process.
Thus, what you ultimately need is to specify a data generating process that you believe is appropriate for the eventual likert data that you want to simulate. After that, there is just the trivial implementational details specific to the software you intend to use (in your case, R). Because people conceptualize likert data as manifest data derived from a latent variable, the most common approach would be to simulate the latent variable according to the theorized distribution (perhaps a normal distribution), and then have a function that maps it to a small ordered set of numbers (e.g., $1, \ldots, 5$). Note that moving from the latent to the manifest variable makes many of the parameters of the latent variable's distribution unidentifiable, so you often needn't bother worrying about them. A simple approach would be to have just the two steps move directly to the final rating, but a more comprehensive approach could model each item with their own set of the two steps, and then have the likert scale combined from the items just as they would be in a real case.
Here is an example, coded in R. I will imagine that there are 5 items that measure the same construct. As such, they are moderately correlated. Two items might be 'reverse scored', but I will assume that doesn't affect the result appreciably so I won't simulate that. However, I will make some more strongly related to the underlying variable than others, and I will make some biased towards higher or lower ratings.
set.seed(8649) # this makes the example exactly reproducible
N = 10 # this is how much data I'll generate
latent = rnorm(N) # this is the actual latent variable I want to be measureing
##### generate latent responses to items
item1 = latent + rnorm(N, mean=0, sd=0.2) # the strongest correlate
item2 = latent + rnorm(N, mean=0, sd=0.3)
item3 = latent + rnorm(N, mean=0, sd=0.5)
item4 = latent + rnorm(N, mean=0, sd=1.0)
item5 = latent + rnorm(N, mean=0, sd=1.2) # the weakest
##### convert latent responses to ordered categories
item1 = findInterval(item1, vec=c(-Inf,-2.5,-1, 1,2.5,Inf)) # fairly unbiased
item2 = findInterval(item2, vec=c(-Inf,-2.5,-1, 1,2.5,Inf))
item3 = findInterval(item3, vec=c(-Inf,-3, -2, 2,3, Inf)) # middle values typical
item4 = findInterval(item4, vec=c(-Inf,-3, -2, 2,3, Inf))
item5 = findInterval(item5, vec=c(-Inf,-3.5,-3,-1,0.5,Inf)) # high ratings typical
##### combined into final scale
manifest = round(rowMeans(cbind(item1, item2, item3, item4, item5)), 1)
manifest
# [1] 3.4 3.6 3.4 3.8 2.6 3.4 3.2 2.0 3.8 3.2
round(latent, 1)
# [1] 1.3 0.6 0.2 1.0 -1.5 0.1 0.4 -2.5 2.3 -0.3
cor(manifest, latent)
# [1] 0.9280074