Rigorous multivariate differentiation of integral with moving boundaries (Leibniz integral rule)

This is probably just a more streamlined version of the proof in the article you linked to, mentioning somewhat less differential forms. I don't know a more elementary proof, but hope this is of use anyway.

I'll write $D_t subset mathbb{R}^d$ for the domain of integration and I'll assume that the map $tmapsto D_t$ is sufficiently regular in the sense that for sufficiently short times, the motion of the domain can be described by the flow of an ambient vector field $X$ on $mathbb{R}^d$. I.e. for each $tin mathbb{R}$ there should exist a vector field $X$ on $mathbb{R}^d$ with flow $Phi_h:mathbb{R}^d to mathbb{R}^d$ such that $$ D_{t+h}=Phi_{h}(D_{t}) qquad forall h in (-epsilon,epsilon). $$

(The field $X$ is not unique, one may add any vector field that's tangent to the boundary of $D_t$.)

With such a vector field $X$ given we compute begin{align} frac{d}{dt} int_{D_t}F(mathbf{x},t)dmathbf{x} &= lim_{hto 0}frac{1}{h}left( int_{D_{t+h}}F(mathbf{x},t+h)dmathbf{x} - int_{D_{t}}F(mathbf{x},t)dmathbf{x}right)\ &=lim_{hto 0}frac{1}{h}left( int_{Phi_h(D_{t})}F(mathbf{x},t+h)dmathbf{x} - int_{D_{t}}F(mathbf{x},t)dmathbf{x}right)\ &=lim_{hto 0}frac{1}{h}left( int_{D_{t}}Phi^*_h(F(mathbf{x},t+h)dmathbf{x}) - int_{D_{t}}F(mathbf{x},t)dmathbf{x}right)\ end{align}

(The previous equality is just the substitution rule for integrals in several variables. So $Phi^*_hleft(F(mathbf{x},t+h)dmathbf{x}right)$ denotes the pullback of the differential form $F(mathbf{x},t+h)dmathbf{x}$ along $Phi_h$. It's what you get by replacing every $mathbf{x}$ with $Phi_h(mathbf{x})$ in the expression $F(mathbf{x},t+h)dmathbf{x}$. We continue by swapping limit and integral)

begin{align} &=int_{D_{t}}lim_{hto 0} frac{Phi^*_h(F(mathbf{x},t+h)dmathbf{x}) - F(mathbf{x},t)dmathbf{x}}{h} end{align}

(To go on use Taylor $F(mathbf{x},t+h)=F(mathbf{x},t)+frac{partial F(mathbf{x},t)}{partial t}cdot h+O(h^2)$ as well as linearity of $Phi^*$ and $Phi^*_0=text{id}_{mathbb{R}^n}$.) $$ =int_{D_{t}} left(frac{partial F(mathbf{x},t)}{partial t}dmathbf{x} +lim_{hto 0}frac{Phi^*_h(F(mathbf{x},t)dmathbf{x}) - F(mathbf{x},t)dmathbf{x}}{h}right). $$ Now that limit in the integral is per definition $L_X(F(mathbf{x},t)dmathbf{x})$, the Lie derivative of the form $omega=F(mathbf{x},t)dmathbf{x}$ with respect to $X$. Using Cartan's magic formula $L_x omega = d(i_X omega)+i_X(domega)$, which in our case simplifies to $L_x omega = d(i_X omega)$ since $F(mathbf{x},t)dmathbf{x}$ is of top degree, we get $$ =int_{D_{t}}frac{partial F(mathbf{x},t)}{partial t}dmathbf{x} +int_{D_{t}} d left(i_X(F(mathbf{x},t)dmathbf{x})right) $$ Applying Stokes theorem to the second integral we finally arrive at $$ frac{d}{dt} int_{D_t}F(mathbf{x},t)dmathbf{x} = int_{D_{t}}frac{partial F(mathbf{x},t)}{partial t}dmathbf{x} +int_{partial D_{t}} i_X(F(mathbf{x},t)dmathbf{x}). $$ For those who want to get rid of that insertion in the second integral, a direct computation shows that $int_Gamma i_X( f(mathbf{x})dmathbf{x})=int_Gamma f(mathbf{x}) Xcdot dvec{A}$ where $dvec{A}$ denotes the usual "surface element" with respect to the standard metric. So you can write the Leibniz integral rule as $$ frac{d}{dt} int_{D_t}F(mathbf{x},t)dmathbf{x} = int_{D_{t}}frac{partial F(mathbf{x},t)}{partial t}dmathbf{x} +int_{partial D_{t}} F(mathbf{x},t), Xcdot dvec{A}. $$