I've just read the other questions and what I have here is basically the full mathematical treatment of what Randy has already said, just to let you know before you plow on!
Standard cartesian rectangular coordinate produce unit vectors [e_x] at each point P in the direction of the positive x-axis and [e_y] at each point P in the direction of the positive y-axis. We may express a position vector at a point P in this basis (the basis of unit vectors [e_x] and [e_y]) as:
[r] = x[e_x] + y[e_y] (1)
If we express vectors in terms of polar coordinates then we first need to express the coordinates of point P (x,y) in terms of polar coordinate (r,θ):
x = rcos(θ) (2i)
y = rsin(θ) (2ii)
Now call [e_r] the radial unit vector pointing away from the origin at each point P and call [e_θ] the angular unit vector pointing (clockwise) tangential to a circle (whose center is the origin) at each point P. We may express the radial position vector as (using (1) and (2)):
[r] = rcos(θ)[e_x] + rsin(θ)[e_y] (3)
The unit radial vector may be expressed by taking the partial derivative of [r] (the position vector) with respect to a change in radial location (whilst holding θ constant) and then dividing by the magnitude of this so that it has a magnitude of 1:
[e_r] = (∂[r]/∂r) / |∂[r]/∂r|
∂[r]/∂r = ∂(rcos(θ)[e_x] + rsin(θ)[e_y])/∂r
=> ∂[r]/∂r = cos(θ)[e_x] + sin(θ)[e_y]
=> |∂[r]/∂r| = √( (∂[r]/∂r).(∂[r]/∂r) )= √(cos^2(θ) + sin^2(θ)) = 1
=> [e_r] = cos(θ)[e_x] + sin(θ)[e_y] (4)
And so we recover something that we already knew:
[r] = r*[e_r] (5)
That is that the position vector at P [r] is given by the radial distance multiplied by the unit radial vector.
We may, however also find the unit angular vector [e_θ] by taking the partial derivative of the position vector [r] (using (3)) with respect to angular displacement whilst holding radius constant, then dividing out it's magnitude as follows:
[e_θ] = (∂[r]/∂θ) / |∂[r]/∂θ|
=> ∂[r]/∂θ = ∂(rcos(θ)[e_x] + rsin(θ)[e_y])/∂θ = -rsin(θ)[e_x] + rcos(θ)[e_y]
=> |∂[r]/∂θ| = √( (∂[r]/∂θ).(∂[r]/∂θ) ) = √(r^2(sin^2(θ) + cos^2(θ))) = r
=> [e_θ] = -sin(θ)[e_x] + cos(θ)[e_y] (6)
It may seem obvious that [e_r] and [e_θ] are perpendicular but we may prove it with the help of (4) and (6):
[e_r].[e_θ] = (cos(θ)[e_x] + sin(θ)[e_y]).(-sin(θ)[e_x] + cos(θ)[e_y])
=> [e_r].[e_θ] = -cos(θ)sin(θ) + sin(θ)cos(θ)
=> [e_r].[e_θ] = 0 (7)
Now we can express the velocity [v] = d[r]/dt in polar coordinates:
[v] = d[r]/dt = d(r*[e_r]) = (dr/dt)*[e_r] + r(d[e_r]/dt) (8)
From (4). Now we may also find d[e_r]/dt:
d[e_r]/dt = (∂[e_r]/∂r)(dr/dt) + (∂[e_r]/∂θ)(dθ/dt) (9)
Using (4) and (6). But of course from (4):
∂[e_r]/∂r = 0 (10)
and
∂[e_r]/∂θ = -sin(θ)[e_x] + cos(θ)[e_y] = [e_θ] (11)
call the angular speed ω = dθ/dt and use (10) and (11) so that (9) becomes:
d[e_r]/dt = ω[e_θ] (12)
Now put (12) into (8):
[v] = (dr/dt)[e_r] + rω[e_θ] (13)
So (13) is our expression for the velocity at a point P using the polar basis vectors [e_r] and [e_θ]. (5) is the expression for the position vector of point P using the same polar basis. Therefore we may take the dot product [r].[v] using (5) and (13) to give:
[r].[v] = (r[e_r]).((dr/dt)[e_r] + rω[e_θ])
=>[r].[v] = r(dr/dt)([e_r].[e_r]) + r^2ω([e_r].[e_θ])
And from (7) along with the fact that [e_r].[e_r] = |[e_r]|^2 = 1 (which is obvious from (4)) we have:
[r].[v] = r(dr/dt) (14)
As you say in the question. Notice that I have mentioned nothing of planets, I have simply expressed position vector [r] and velocity vector [v] at a point P(x,y) or (r,θ) in polar coordinates. So (14) is true regardless of the scenario. As the other answers here say, the orbit of planet is ALWAYS elliptical. In high school physics we tend to approximate to a circle but there is always some orbit eccentricity no matter how small. Thus dr/dt ≠ 0 for planetary orbits ever in reality, only in simplified models.