Question on using Leibniz formula to derive thin-film equation from Navier-Stokes I actually posted this to math.stackexchange.com a few months ago but never got any answers.
I am trying to work through the derivation in this paper by Petr Vita, which derives a thin-film simplification of the Navier-Stokes equation, similar to the Reynolds or Lubrication Equation, but including inertial terms as well.  To walk through the major steps to the point where I have questions:


*

*Start with basic N-S equation (eq 1 in the paper):
$$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot(\mathbf{u u}) \right) = -\nabla p + \rho \mathbf{g} + \nabla \cdot \mathbf{\underline{T}} $$
Here $ \mathbf{\underline{T}}$ is the deviatoric stress tensor, and I have left off the final body force term since it isn't used in the rest of the paper.

*Use the thin-film assumption that $u_z=0$ and define (eq 4 in the paper) $$\bar{\mathbf{u}}=\frac{1}{h}\int_0^h \mathbf{u}\, dz$$

*Integrate the N-S equation with respect to $z$ from $0$ to $h$ (equation 8 in the paper). $$ \rho \frac{\partial}{\partial t} (h \bar{\mathbf{u}}) + \rho \int_0^h \nabla \cdot (\mathbf{u u})dz = -h \nabla p -\left.\mu \frac{\partial \mathbf{u}}{\partial z}\right|_{z=0} $$
Obviously there are a few steps being skipped over here.  For the time-derivative term, I use the Leibniz formula to derive the following:
$$ \int_0^h \rho \frac{\partial \mathbf u}{\partial t} dz = \rho \frac{\partial}{\partial t} \int_0^h \mathbf u\,dz - \rho \frac{\partial h}{\partial t} \mathbf u (x,y,h,t) + \rho \frac{\partial 0}{\partial t} \mathbf u (x,y,0,t)$$
Obviously the last term is $0$ and can be dropped. Also $ \rho \frac{\partial}{\partial t} \int_0^h \mathbf u\,dz = \rho \frac{\partial}{\partial t} (h \bar{\mathbf{u}}) $, giving the form seen in the equation.  However, I don't see how $ \color{blue}{\rho \frac{\partial h}{\partial t} \mathbf u (x,y,h,t)} $ can be taken to be $0$.  The height of the film is certainly changing with time, and the top surface has a von Neumann boundary condition, not a Dirichlet no-slip boundary. Any insight here?
Also, the deviatoric stress has to be integrated as well.  I think the divergence theorem can be used here:
$$ \int_V \nabla \cdot \mathbf{\underline{T}}\, dV = \int_S \mathbf{n} \cdot \mathbf{\underline{T}} \, dS$$
In this case that should come out to be
$$ \int_0^h \nabla \cdot \mathbf{\underline{T}}\, dz =\left.\mu \frac{\partial \mathbf{u}}{\partial z}\right|_{z=h} - \left.\mu \frac{\partial \mathbf{u}}{\partial z}\right|_{z=0} $$
The top surface stress is $0$, leaving the bottom stress term as is found in the derived equation, right?

*Now we get to my main question, the integration of the $\nabla \cdot(\mathbf{u u})$ term.  The author is able to evaluate this term by using a combination of Pohlhausen's method of assuming a cubic profile for the liquid flow, and the Reynold's Averaged Navier-Stokes method of splitting the velocity into an average velocity and deviation from that average.
For the cubic profile he defines:
$$ \mathbf{u}(x,y,z) = u(x,y,\xi) \text{, where} $$
$$ u(x,y,\xi) = a_0 + a_1\xi + a_2\xi^2 + a_3\xi^3,\quad \xi \in \langle 0,1
\rangle,\; z=h\xi $$
Then he applies the boundary conditions and integral relation to obtain 
$$ u(x,y,\xi) = \mathbf{u}_{disk} + (\bar{\mathbf{u}}-\mathbf{u}_{disk})\left( \frac{12}{5}\xi - \frac{4}{5}\xi^3 \right) $$
This step is fine, I had no problems figuring it out.  Then the author defines the velocity fluctuation (with respect to the vertical direction) $\mathbf{\tilde u}$ as 
$$ \mathbf{u} = \mathbf{\bar u} + \mathbf{\tilde u} \text{.  This makes:} $$
$$ \int_0^h \mathbf{\bar u}\, dz = h\,\mathbf{\bar u}\text{, and } \int_0^h \mathbf{\tilde u}\, dz = 0 \text{.}$$
Anyway, so the author does this to integrate the advection term:
$$ \int_0^h \nabla \cdot ( \mathbf{uu} )\,dz = \nabla \cdot \left( \int_0^h \left[ \mathbf{\bar u} \mathbf{\bar u}+\mathbf{\bar u} \mathbf{\tilde u} + \mathbf{\tilde u} \mathbf{\bar u} + \mathbf{\tilde u} \mathbf{\tilde u} \right]\,dz \right)$$
So on the RHS the 1st term is just a constant, the 2nd and 3rd terms become $0$, and then he uses the derived polynomial form of $u$ to evaluate the last term.  However, how was he able to pull the divergence operator out of the integral?  He didn't use the divergence theorem, and I don't know if you can use the Leibniz formula on a divergence operator.  If you could do that though, wouldn't you have a term that's something like $\nabla \cdot h \left.(\mathbf{u}\mathbf{u})\right|_{z=h}$ leftover as well? 
 A: I think that most problems come from the fact, that you do not consider film height $h$ as a function $h(x,y,t)$ whole time:
$$
\int_{a(\theta)}^{b(\theta)}\partial_\theta f(\theta, z) dz = \frac{d}{d\theta} \left(\int_{a(\theta)}^{b(\theta)} f(\theta, z) dz\right) + f(\theta, a(\theta))a'(\theta) - f(\theta, b(\theta))b'(\theta)
$$
The integral over film depth does not collide with film height function parameters.
I would be careful with using term Reynold's Averaged Navier-Stokes (RANS) in this case. The splitting of a quantity into average and its deviation is Reynold's Decomposition. RANS is a product of Reynold's Decomposition applied to velocity in time, where here is the decomposition done in space. Warning: Term RANS makes many members of fluid-dynamics community blind to your idea straight-away.
I used the word 'quantity' for a reason. You can use the technique from the paper for any quantity, really (with exception of pressure, I think). This can be very useful if you want to bring energy or chemistry model into your system.
A: Recall the objection to taking the derivative operator outside an integral in the first place - you can't do it when the domain of integration is itself a function of what you're trying to differentiate with respect to, e.g. time.
In this paper, however, $\int_{h}$ seems to be an integral over a spatial region that is only a function of time, namely the thin film itself. So $h$ is not a function of space (x,y,z) (if I read it right).
So  you can think of this as a generalized Leibniz rule application (assuming sufficiently continuous derivatives etc), i.e. a generalization of
$$
\int_{a}^{b} \frac{\partial f(z,x)}{\partial x} dz = \frac{d}{dx} \int_{a}^{b} f(z,x) dz
$$
when $a,b$ are independent of $x$. If $a, b$ depended on $x$ you'd have two other terms in the Leibniz equation. 
PS: Note that normally you use the Leibniz rule and its generalizations to evaluate the derivative of an integral, but there is nothing that prevents you from going the other way, which is why I switched the left and right hand sides in writing it.
PS2: Oh wait, further down in the paper they have h as a function of x and y as well. Bummer. Let me think about this and get back to you. 
