The continuous model of the problem suggests that the limit does depend on $f$ (and $u$). More precisely, it depends on how fast the parameter $f$ is suppressed in the expression whose limit you are taking;  behavior of $\lim y(nt, f_n)/n$ will depend on the limit of $nf_n$ as $n \to \infty$. The answer will be a function of the limit of $nuf_n$. Only when this limit is zero does one get the proposed formula.

The associated differential equation is $Y' = 1/(1+Ae^{BY})$ where $A = e^z$ and $B=uf$. Its solution vanishing at 0 is $Y(t) = H^{-1}(t)$ where $H(t) = t + (A/B)(e^{Bt} - 1)$.  It does look like this matches the asymptotic behavior of your sequence for $y$ both in the large and small range.  For small $t$, the expansion $Y(t) = t/(1+A) + O(t^2)$ corresponds to your formula, but I think the answer is not quite that simple: you have to establish whether the expression whose limit is calculated belongs to the small regime where $Y(t)$ is approximately linear, or the large regime where $Y(t)$ is logarithmic, $Y(t)=O(\log(t))$.  The limit uses $n$ iterations so we want to know, as a function of $B \sim 1/n$, whether the transition between regimes happens at a point much larger than $1/B$.   However, it's easy to calculate that the ratio $H(t)/t$ moves away from 1 (the difference is larger than some constant independent of $B$) as soon as $Bt$ is of order 1 (i.e., bounded below by a given positive constant) and this would spoil the limit if the differential equation is a good model of the difference equation. 

(ADDED: for comparison of $Y$ predictions with $y$ simulations, in the phase transition where $Bx$ is of order 1, $Y(x) \sim x/C$ and $H(x) \sim Cx$, with $C = 1 + A(e^q - 1)/q$,  and $q = Bt = ufx$.  That is, $Y$ stays approximately linear but the coefficient goes to zero, consistent with the idea that it's turning into a logarithmic function.   In the notation of the question, $x=nt, \quad f=f_0/n$, for some constant $f_0$ and with $u$ and $t$ also held constant while $f$ varies with $n$, so that the phase transition parameter is $q=uf_0 t$ and the predicted value of the limit, if $Y$ is a good approximation for $y$, is $L_{pred} = \lim Y(nt)/n = (q/(q + Ae^q - A))$.  Does this match the simulations?)

To see the small-$y$ behavior directly in the difference equation, it can be expanded in powers of $y$. 

$y(t+1) - y(t) = 1/(1+A)  - (AB)/(1+A)^2)y + O(y^2)$

Your formula proposes that when $B \sim 1/n$, the effect of the $y^{\geq 1}$ terms is of order smaller than $n$ for $t \in [0,n]$.   The sum of the first $t$ values of the $y^1$ term will be of order $t^2$, so one expects these corrections to be suppressed only on a short interval, $t << n^{1/2}$.   The calculation with the differential equation suggests that $f_n = f/n$ is too large a parameter ; this calculation with the truncated difference equation can be used to prove that $f_n = f/n^k$ is small enough for any $k > 2$.  Adding higher degree terms to the approximate difference equation would, I suppose, only get closer to the picture suggested by the differential equation. 

To prove rigorously the predictions from the differential equation you could try to control $y$ by trapping the sequence $y(n)$ between two trajectories of the ODE.  If simulations are consistent with a heuristically "wrong" formula it would be very interesting to sort out what the truth is.