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 $Bt$ is of order 1, $Y(t) \sim t/C$ and $H(t) \sim Ct$, with $C = 1 + A(e^q - 1)/q$, and $q = Bt =uft$. That is, $Y$ stays approximately linear but the coefficient goes to zero, consistent with the idea that it's turning into a logarithmic function. Let $t=nt_0, \quad f=f_0/n$, for some constant $f_0$ and with $u$ and $t_0$ also held constant while $f$ varies with $n$, so that the phase transition parameter is $q=uf_0 t_0$ and the predicted value of the limit, if $Y$ is a good approximation for $y$, is $L_{pred} = \lim Y(nt_0,f_0/n)/n = \lim nt_0/nC = t_0/C = t_0(q/(q + Ae^q - A))$. In the original notation of the question, $L = t/(1+{e^z}F(uft))$ where $F(x)=(e^x-1)/x$. 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.
