But $f(N)$ can be any complicated function. How can we show: $𝑆\sim𝑁$?
The answer is: you can't, because it's not, in general. Even in your own expression, with
\begin{align}
\Omega = f(N) V^N U^{Nf/2}
\end{align}
you have to ask: how do $V$ and $U$ scale with $N$? That's what the people pointing to "intensive" versus "extensive" variables are talking about. For example, if I double $N$ by copying the system and setting the copy next to it, how does the variable change (likewise, I could cut it in half)?
Tl;dr: you have to assume that the particles that make up the system are weakly interacting/independent. If you can do that, then you can derive the constraints on the form of $f(N)$, and those constraints will lead to $S\sim N$.
Not a full derivation, but here are some of the specifics:
The definition of entropy was given by Gibbs:
\begin{align}
S = -k_B \sum_{i} P_i\ln P_i,
\end{align}
where $P_i$ is the probability that the system is in the state labeled by $i$. Let's split the system into two parts, and give the two parts the labels $j$ and $l$. The formula becomes:
\begin{align}
S = -k_B \sum_{jl} P_{jl}\ln P_{jl},
\end{align}
where $P_{jl}$ is what's known as a joint probability: the probability that the first state label that we can call $J$ has the value $j$, and the second state label $L$ has the value $l$ at the same time. Note that $\sum_{jl} P_{jl} = 1$, because probabilities must sum to 1, by definition.
There's something called conditional probability $P_{j|l}$ where we read this as, "The probability that $J=j$ (randomly) given that $L=l$ (non-randomly)." In equations it is defined by the relation to the joint probability by
\begin{align}
P_{jl} = P_{j|l} P_l,
\end{align}
where $P_l$ is called the marginal probability and it is defined by $P_l = \sum_j P_{jl}$.
Last definition. If $P_{j|l} = P_j$ (i.e. the conditional probability for $j$ doesn't depend on $l$) then we say that the random variables $J$ and $L$ are independent.
Now, let's put our definitions to use. First, we insert the definition for conditional probability into the formula for entropy:
\begin{align}
S &= -k_b \sum_{jl} P_{j|l}P_l \ln\left(P_{j|l}P_l\right) \\
&= -k_b \sum_{l}P_l\left(\sum_j P_{j|l}\left[ \ln P_{j|l} + \ln P_l\right]\right) \\
&= -k_b \left[\sum_{l}\left(P_l\sum_j P_{j|l} \ln P_{j|l}\right) + \sum_{l}\left(P_l\ln P_l\sum_j P_{j|l}\right)\right] \\
&\qquad \text{because }\sum_{j} P_{j|l} = 1: \\
&= -k_b \left[\sum_{l}\left(P_l\sum_j P_{j|l} \ln P_{j|l}\right) + \sum_{l}P_l\ln P_l\right] \\
&= S(j|l) + S(l),
\end{align}
where the last line is the definition of the conditional entropy $S(j|l)$. Note that if $P_{j|l} = P_j$ then $S(j|l) = S(j)$.
The necessary and sufficient condition for $S\sim N$ is this: if we label the states of the $N^{\mathrm{th}}$ particle with $l$ and the states of every other particle with $j$, then $P_{j|l} = P_j$, and that this holds for every $N$ between $1$ and the $N$ of interest (proof by induction). In other words, $S\sim N$ (entropy is extensive in $N$) when you can assume that the interactions between the particles are weak enough to be negligible.
It's an assumption that works very well most of the time. It's what we assume for the ideal gas. For solids, where the molecules interact quite strongly, you can either talk about how the number of particle on the surface boundary between divisions of the system in half are much fewer than the number of molecules in the bulk, or you can say that the particles of interest are the phonons (sound equivalent to photons) that travel through the solid as non-interacting waves.
The assumption runs into trouble when you can't make these kind of approximations. For example, a Coulomb gas (i.e. a plasma) has long-range interactions that make it a non-ideal gas and its entropy may not be extensive. You can also add long-range interactions to the "simple" Ising model. The results of doing that are being actively studied and are known to produce "entanglement entropy," the quantum mechanical version of non-independent conditional probability.
