I'm numerically solving a system of PDE's consiting of some conservation laws (Euler equations) along with advection of non-additive values (like molar mass or unit heat capacity, for example). The question is concerned with advection equations. In 2D cartesian coords, advection of value $\phi(x,y,t)$ is governed by this equation:
$$
\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} + v \frac{\partial \phi}{\partial y} = 0
$$
where $u$ and $v$ are velocity components along $x$ and $y$ axes. The velocity vector is provided by solving Euler equations.
This equation may be written in quasi-conservative form with source term:
$$
\frac{\partial \phi}{\partial t} + \frac{\partial \phi u}{\partial x} + \frac{\partial \phi v}{\partial y} = \phi\frac{\partial u}{\partial x} + \phi\frac{\partial v}{\partial y}
$$
In this form it is being successfully solved by some numerical method. By "succesfully" I mean that the solution looks physical enough: values of $\phi$ are advected with the gas flow and are not "accumulated" in the areas where gas flow converges (as would gas density accumulate, for example).
Now I try to solve the same system in axysimmetric cylindrical coordinates $(x, r)$. I found this form of 1D scalar advection equation in the literature:
$$
\frac{\partial \phi}{\partial t} + \frac{v}{r} \frac{\partial r\phi}{\partial r} = 0
$$
Or, in quasi-conservative form with source term (also, multiplied by $r$):
$$
\frac{\partial r\phi}{\partial t} + \frac{\partial r \phi v}{\partial r} = r \phi\frac{\partial v}{\partial r}
$$
(I derived this question by hand, so cannot 100% guarantee that is correct, though)
One exact solution is given in the paper:
$$
\phi(r,0) = \left\{\begin{array}{ll}
\frac{sin^4(\pi r)}{r}, & 0\leq r\leq 1\\
0, & r > 1
\end{array}\right\},
$$
$$
\phi(r,t) = \left\{\begin{array}{ll}
\frac{sin^4(\pi (r-vt))}{r}, & vt \leq r\leq vt + 1\\
0, & \text{otherwise.}
\end{array}\right\}
$$
Here is graphic representation of this solution for $v = 1$ at $t = 1$ from that paper:
It seems that original function shape is moved by 1 unit of distance away from symmetry axis and diminished in height by a factor of about ~3.
My questions are:
- Why does the height of the function shape diminish? In Cartesian coords, the original shape would just move and preserve it's shape. To my intuition, such phenomenon looks more like a conservation law behaviour in polar coords: function height decreases due to its base expansion so that total integral of the function remains constant.
- How should one properly describe an advection (or transport) of some non-additive value in cylindrical coords? I.e. of such a value that is simply "glued" to fluid particles and would not decrease in such divergent flow as shown above (or increase in convergent flow).
Best Answer
Yes, the above observations and derivations seem correct.
In cylindrical coordinates, the advection equation $\partial_t \phi + \frac{v}{r} \partial_r (r\phi) = 0$ for a given quantity $\phi(r,t)$ transported with uniform velocity $v\equiv 1$ m/s along the radial direction is a conservation law. In fact, it might be rewritten as $$ \partial_t \phi + \text{div}(\phi \boldsymbol v) = 0 \, , \qquad \boldsymbol v= v \hat{\boldsymbol r} $$ where $\hat{\boldsymbol r}$ is the radial unit vector. Now, integration yields $$ 0 = \int_{\Omega} \left[\partial_t \phi + \text{div}(\phi \boldsymbol v)\right] \text{d}V = \frac{\text d}{\text d t} \int_{\Omega} \phi\, \text{d}V + \int_{\partial\Omega} \phi \boldsymbol v \cdot \boldsymbol{n}\, \text{d}S \, , $$ where we have used the divergence theorem. Here, $\boldsymbol n$ is the outgoing normal unit vector to the boundary $\partial\Omega$ of the spatial domain $\Omega$. Therefore, the variation in time of the total quantity $\int \phi \, \text{d}V$ can be expressed in terms of the flux $\phi \boldsymbol v$ passing through the boundary.
If in-plane transport along a given direction is desired, a straightforward option would be to transform the initial Cartesian transport equation (a.k.a. the "color equation") $$ \partial_t \phi + \text{div}(\phi \boldsymbol{v}) = 0, \qquad \boldsymbol{v} = u\hat{\bf x}+v\hat{\bf y} $$ with $u$, $v$ constant, into cylindrical coordinates. This is done by using the relationships \begin{aligned} \hat{\bf x} &= \cos\theta\, \hat{\boldsymbol r} - \sin\theta\, \hat{\boldsymbol\theta} \\ \hat{\bf y} &= \sin\theta\, \hat{\boldsymbol r} + \cos\theta\, \hat{\boldsymbol\theta} \end{aligned} between unit vectors. Therefore, the velocity field becomes velocity $\boldsymbol{v} = v_r \hat{\boldsymbol r} + v_\theta \hat{\boldsymbol \theta}$ , and the transport equation $$ \partial_t \phi + v_r \partial_r \phi + \frac{v_\theta}{r} \partial_\theta \phi = 0 $$ $$ \text{with}\qquad v_r = u\cos\theta + v\sin\theta ,\qquad v_\theta = -u\sin\theta + v\cos\theta $$ follows from the expression of the divergence operator in cylindrical coordinates.
Note: The product rule for the divergence operator reads (cf. vector calculus identities) $$ \text{div}(\phi \boldsymbol{v}) = \phi \,\text{div}\, \boldsymbol{v} + (\text{grad}\, \phi)\cdot \boldsymbol{v} \, . $$ By definition, the material derivative $\dot\phi$ in Eulerian coordinates reads $$ \dot\phi = \partial_t \phi + (\text{grad}\, \phi)\cdot \boldsymbol{v} $$ so that the conservation law $\partial_t \phi =- \text{div}(\phi \boldsymbol{v})$ amounts to the evolution equation $$ \dot \phi/\phi = -\text{div}\, \boldsymbol{v} \, . $$ Therefore, the quantity $\dot \phi$ vanishes in incompressible or in uniform flow, in which cases the advection equation $\dot\phi = 0$ is equivalent to the conservation law. In particular, note that the conservation law is satisfied by the mass density $\phi = \rho$, owing to the balance principles of continuum physics (cf. integral form above).