Momentum of slowly spinning (viscous) fluid If we have a massless cylindrical container (or radius $R$) with a liquid of certain density $\rho$ and viscosity $\mu$ at rest. Then at time zero we impart a constant rotational velocity $\Omega$ on the cylinder and watch as the liquid accelerates from the outer walls inwards (due to the viscosity).
I want to know what is the function of total angular momentum $L$ with time (and consequently the effective mass moment of inertia of the liquid $L=I_{eff} \Omega$). More specifically the radius of gyration $I_{eff} = m \kappa^2$ where $m=\rho \pi R^2 h$.
I looked at concentric cylindrical slices in order to derive the equations of motion but I am stumbling at the shape of the tangential velocity as a function of radius $r$ and time $t$.
Based on of $\nabla$ in cylindrical coordinates (for a Newtonian fluid) I think shear stress is $$\tau_{r\theta} = \mu \left( \frac{\partial v_\theta}{\partial r} - \frac{v_\theta}{r} \right)$$ 
The cylindrical slice has surface area $A = 2 \pi r h$, or ${\rm d}A = 2 \pi h {\rm d}r$. 
The volume is ${\rm d}V = A {\rm d}r  = 2 \pi r h {\rm d}r$. 
I think the radial force balance acting on the volume due to the shear stress is
$$ A {\rm d}\tau_{r\theta} + \tau_{r\theta} {\rm d}A = \dot{v_\theta} \rho {\rm d} V$$
In which this leads to a differential equation $$\tau_{r\theta}'= \frac{\partial \tau_{r\theta}}{\partial r} = \rho \dot{v_\theta} - \frac{\tau_{r\theta}}{r} $$
The shear stress slope is (from chain rule) $$ \tau_{r\theta}' = \mu \left( v''_{\theta} - \frac{v'_\theta}{r}+ \frac{v_\theta}{r^2} \right) $$
With some algebra I get finally that
$$\dot{v}_{ \theta} = \frac{\mu}{\rho} v''_{\theta} $$
So the acceleration of the cylindrical slice is proportional to the curvature of the velocity profile.
Here is where I am stuck. I am not sure how to proceed to derive $v_\theta (r,t)$.
 A: Assuming the equation you found is correct, it's just the 1-dimension heat equation:
$$
D^2 \frac{\partial^2 v}{\partial r^2} = \frac{\partial v}{\partial t},
$$
where $D^2 \equiv \mu/\rho$.  We want to solve it on the domain $r \in [0, R]$, $t \in [0, \infty)$, subject to the boundary conditions $v(r,0) = 0$, $v(R, t) = R \Omega$, and $v(0,t) = 0$.  (This last condition is inserted so that the velocity field remains continuous at the origin;  remember that $v$ is only the tangential component of the field.)  
The steady-state solution $v_\infty(r,t)$ to this equation is pretty obvious:  if $\dot{v}_\infty = 0$, then $v''_\infty = 0$ as well, and so
$$
v_\infty(r,t) = \Omega r.
$$
This makes sense:  the whole cylinder is rotating rigidly at late times.
Now define $\delta v(r,t) \equiv v(r,t) - v_\infty(r,t)$.  By construction, $\delta v$ also satisfies the heat equation, but with different boundary conditions:
$$
\delta v(R,t) = \delta v (0, t) = 0;  \qquad \delta v (r,0) = - \Omega r \equiv \delta v_0(r).
$$
This is just a standard heat-diffusion problem with an initial heat distribution $\delta v_0(r)$.  To solve this, we use separation of variables;  it's not too hard to show that any solution of the form
$$
f_n(r,t) = \sin \left( \frac{n \pi r}{R} \right) e^{-D^2 \pi^2 n^2 t/R^2}
$$
will work.  Assuming that $\delta v(r,t) = \sum_n A_n f_n(r,t)$, we have
$$
\delta v_0 (r) = - \Omega r = \sum_n A_n \sin \left( \frac{n \pi r}{R} \right),
$$
i.e., $\delta v_0$ is expressed as a Fourier series. Working through this (I used Mathematica to save time), we get
$$
A_n = \frac{2 \Omega R (-1)^n}{\pi n}
$$
Thus, the full solution for the tangential velocity as a function of time and space is:
$$
\boxed{ v(r,t) = \Omega r + 2 \Omega R \sum_{n=1}^\infty \frac{(-1)^n}{\pi n} \sin \left( \frac{n \pi r}{R} \right) e^{-D^2 \pi^2 n^2 t/R^2}.}
$$
It may be possible to sum this up into a closed-form expression;  I'll let you know if I make any progress on this.  
A: First of all, I think the equation of motion is not correct. I believe the equation of motion is 
$$
\dot{v}_\theta  = \nu \left( v^{\prime\prime}_\theta + \frac{v^\prime_\theta}{r} - \frac{v_\theta}{r^2} \right) ,
$$
where $\nu=\mu/\rho$. This is easiest to derive using the result for the Laplacian of a vector field in cylindrical coordinates. 
In order to derive this result from force balance we have to consider a volume element $V=L_z L_\theta L_r$ in cylindrical coordinates. Here, $L_\theta=r\Delta\theta$ and $L_r=\Delta r$. We can ignore $L_z=\Delta z$ as nothing depends on $z$. The force on $L_\theta$ is $F=\tau_{r\theta}\hat{e}_\theta L_\theta$. The net acceleration is due to the difference between the forces on the faces at $r$ and $r+\Delta r$. We get $\dot{v}_\theta=\nu[\tau_{r\theta}'+\tau_{r\theta}/r]$ as explained in the question. The force on $L_r$ is $F=\tau_{r\theta}\hat{e}_r$, and the net acceleration is due to the difference between the forces on $L_r(\theta)$ and $L_r(\theta+\Delta\theta)$. Using $\partial \hat{e}_r/(\partial \theta)=\hat{e}_\theta$ we get $\dot{v}_\theta=\nu\tau_{r\theta}/r$. Combining both we get
$$
\dot{v}_\theta=\nu\left(\tau_{r\theta}'+2\frac{\tau_{r\theta}}{r}\right),
$$
which agrees with the formula above.  
The basic time scale is given by vorticity diffusion, so we expect $v_\theta \sim \exp(-c\nu t/R^2)$, but a little more effort is required to get $c$ and the precise $r$ dependence. You can tackle this by Fourier expanions, but because of the cylindrical geometry you should really use a Bessel expansion. 
Make a separation ansatz $v_\theta(r,t)=g(r)f(t)$. The equation of $g(r)$ is the Bessel $J_1$ differential equation. Then we expand the boundary condition in $J_1(\lambda_n r/R)$, using orthogonality of the Bessel function with respect to the zeros $\lambda_n$. Then (for spinning down)
$$
v_\theta = 2\Omega R \sum_n \frac{J_1(\lambda_n r/R)}{\lambda_nJ_0(\lambda_n)}
 \exp\left(-\lambda_n\nu t/R^2\right).
$$
The solution for spinning up is just $\exp()\to 1-\exp()$.
Bonus: To convince myself that the solution does indeed satisfy the boundary conditions, I plotted the answer (oscillations near the boundary are due to slow convergence of the Bessel expansion). The figure shows $v_\theta(r)$ for different $t$.
 
A: Edit: This use of non-dimensionalization is wrong because it is incompatible with the boundary conditions. More information here.
We're looking for a solution to the fluid velocity $\vec u$ that looks like $\vec u = u(r,t)\hat\theta$. Given the symmetries of the problem, we have to solve Navier Stokes:
$$
\partial_t u=\nu(\nabla^2 u - \frac u{r^2})
$$
The only quantities that appear in this equation, apart from $u$, are $t$, $r$ and $\nu$. Because $u$ has to be a function of non-dimensional variables, we look for non-dimensional combinations of those quantities and find only one: $\frac{r^2}{\nu t}$. We'll call this quantity $z$ and simplify our expression for $\vec u$:
$$
\vec u = u(z)\hat\theta
$$
Now let's look at the derivatives. Because $u(z)$ can only depend on $r$ and $t$ through $z$, we have
$$
\partial_t u = \partial_t z u' = -\frac z t u'\\
\partial_r u = \partial_r z u' = 2\frac z r u'
$$
This leads to the following expression for the laplacian:
$$
\nabla^2 = \frac 1 r \partial_r(r\partial_r)\\
= \frac{2z}{r^2}\frac d{dz}\left(2z\frac d{dz}\right)\\
= \frac{4z}{r^2}\left(\frac d{dz}+z\frac{d^2}{dz^2}\right)
$$
Inserting this into Navier-Stokes, we get
$$
-\frac z t u' = \frac{4z\nu}{r^2}(u'+zu'')-\nu\frac u{r^2}
$$
We multiply everything by $\frac{r^2}{4z^2\nu}$:
$$
-\frac 1 4 u' = u''+\frac{u'}z-\frac u{4z^2}
$$
Rearranging,
$$
u''+\left(\frac 1 4+\frac 1 z\right)u'-\frac 1{4z^2} u = 0
$$
That's all I have so far, but you can see that through non-dimensionalization you can turn the original PDE into a (simpler?) ODE.
A: I didn't forget about this one! I cannot find fault with your derivation and the resulting PDE. I may have missed something also though, so who knows. But what I can do is come up with a solution for the PDE you do have, which should help if it turns out the original isn't right for one reason or another.
Starting from (where $\nu = \mu/\rho$):
$$ \frac{\partial v_\theta}{\partial t} = \nu \frac{\partial^2 v_\theta}{\partial r^2} $$
and taking the Fourier transform from $r\rightarrow k$, we get:
$$ \frac{\partial \hat{v}_{\theta}}{\partial t} = \nu (i k)^2 \hat{v}_{\theta}$$
which has an exponential function as a solution. The solution is thus:
$$ \hat{v}(k,t) = e^{-\nu k^2 t} \hat{v}(k,0) $$
where $\hat{v}(k,0)$ is the Fourier transform of your initial conditions (the impulsively started wall at $r = R$). For those initial conditions, assuming the wall is started with velocity $U$, the transform of the initial conditions is:
$$\hat{v}_\theta (k,0) = \frac{1}{\sqrt{2\pi}} e^{i k R} U $$
We can now take the inverse Fourier transform of this result and we will get:
$$ v_\theta (r,t) = \frac{U}{2\sqrt{\pi \nu t}} e^{-(r-R)^2/4\nu t}$$
Or at least that's what the math says. The part that bothers me now are the units; the RHS seems to have units of $\text{s}^{-1}$ unless there are some hidden units in $\pi$, which given the weirdness of Fourier transforms wouldn't surprise me at all. 
A: Question; 

I am not sure how to proceed to derive $vθ(r,t)$

It may be easier to make a clean simple start rather than trying to seek a mistake. 
If the fluid is rotating slowly, and it's viscous. You basically have a laminar flow conditions, which makes the whole really simple as the velocity distribution is an exponential curve. Ie parabel. At your case the whole flow will finally rotate as a rigid object, but up to there the velocity change is a cut from this parabel curve.
Here's a helping video about "Rotating Flows".. The needed stuff is at 1-4 min. Basically this picture; 

$U$ is the velocity. The stuff at 3rd row is the specialties in rotating flows; Centrifugal force and Coriolis force. 
I Hope this short answer helps you to proceed.
