# Question on using Leibniz formula to derive thin-film equation from Navier-Stokes

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:

1. 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.
2. 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$$
3. 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?

4. 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?