Note that $v_1$ and $v_2$ are equal in magnitude and perpendicular to the respective radii. So they form an isosceles triangle and the angle between them is the same as the angle between the radii. The base angles of the two isosceles triangles are also equal, hence they are similar.
To show the angles are equal translate the configuration at $P_2$ to $P_1$ and work round the angles in that figure. The angle between $v_2$ and $OP_1$ adds to the angle between $OP_2$ and $OP_1$ to give a right angle and similarly with the angle between $v_2$ and $v_1$ so the angles you want are equal.
The points on the infinite right circular cylinder are at distance $r$ from the axis of the cylinder. The distance between point $\vec{p}$ and line through $\vec{p}_1$ and $\vec{p}_2$ is
$$d = \frac{\left \lvert (\vec{p} - \vec{p}_1) \times (\vec{p} - \vec{p}_2) \right \rvert}{\left \lvert \vec{p}_1 - \vec{p}_2 \right \rvert}$$
so we can square that to get our equation describing an infinite right circular cylinder of radius $r$ whose axis passes through points $\vec{p}_1$ and $\vec{p}_2$:
$$r^2 = \frac{\left \lvert (\vec{p} - \vec{p}_1) \times (\vec{p} - \vec{p}_2) \right \rvert^2}{\left \lvert \vec{p}_1 - \vec{p}_2 \right \rvert^2}$$
That can be written using Cartesian coordinates as
$$r^2 = \frac{ \left ( (y-y_1)(z-z_2) - (z-z_1)(y-y_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 } + \frac{ \left ( (z-z_1)(x-x_2) - (x-x_1)(z-z_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 } + \frac{ \left ( (x-x_1)(y-y_2) - (y-y_1)(x-x_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 }$$
which we can convert to standard form by multiplying by the right-hand-side denominator, and moving all to the same side:
$$\left ( (y-y_1)(z-z_2) - (z-z_1)(y-y_2) \right )^2 + \left ( (z-z_1)(x-x_2) - (x-x_1)(z-z_2) \right )^2 + \left ( (x-x_1)(y-y_2) - (y-y_1)(x-x_2) \right )^2 - r^2 \left ( (x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 \right ) = 0$$
I have at my Wikipedia talk page the equations to locate the intersection between a line passing through origin and an arbitrary right circular cylinder; without caps, with flat end caps, or with spherical end caps. The key point is that when the ray (or line we are testing intersection for) passes through origin (so we use an unit vector $\hat{n}$, $\lvert\hat{n}\rvert=1$, to describe the line or ray), the equations become rather straightforward.
If the underlying problem is related to grinding spheres using a grinding wheel -- a cutting plane, essentially -- then a simplified coordinate system where origin is the sphere origin, and $\hat{z}$ is along the rotation axis, makes for easy calculations. Consider this illustration:
The contact point scribes a circle around the axis of rotation. We can describe this circle using a single scalar, $d$, which is the signed distance from the origin to the plane of that circle, along the axis of rotation. If the circle is not on the surface of the sphere, but perhaps inside it, then we need a second scalar, $r$, to describe the radius of the circle (around the axis of the rotation).
If the radius of the sphere is $R$, then on the surface of the sphere $r = \sqrt{R^2-d^2}$, of course.
If we include the cutting or abrasion done by the wheel, then two pairs of scalars are needed: $(d_1,r_1)$ and $(d_2,r_2)$. If we assume the grinding wheel is planar, and stays at a fixed location with a perfect cut for at least one rotation of the sphere, the removed section of the sphere leaves a flat facet. Mathematically, $$r(d) = r_1 + \left ( d - d_1 \right ) \frac{r_2 - r_1}{d_2 - d_1}$$
where $d_1 \le d \le d_2$. Remember, $d$ is the distance along the axis of rotation, and $r$ the circle radius (around a point on the axis of rotation, at distance $d$ from the center of the sphere); that is why the dependence above is linear too.
If the axis of the grinding wheel always intersects the rotation axis, and you call up $z$ and right $x$ in the above diagram, then $(d,r) = (z,x)$, and the primary problem (surface of revolution caused by grinding, or shape of the "sphere" after several grinds) reduces to planar (2D) vector calculus.
If we also have the axis of rotation $\hat{n}$ in some fixed sphere coordinate system, it is easy to calculate the "ribbons" each grind cuts into the sphere, either parametrically (using say $0 \le u \le \pi$ along the circular ribbon, and $0 \le v \le 1$ across the ribbon -- useful for visualization) or algebraically.
Best Answer
There is a slight mistake in your analysis towards the beginning. You assumed that the speed is equal to the total acceleration time the time. However the radial acceleration has no effect on the magnitude of the velocity in circular motion.
Therefore we can write,
$$ v(t) = a_t t $$
Now lets assume your target velocity is $v_{max}$ then we know ahead of time that the maximum radial acceleration will be,
$$ a_{r_{max}} = \frac{v_{max}^2}{R}$$
Since we also know the maximum total acceleration, $a_{max}$ we can easily compute the maximum allowed tangential acceleration.
$$(a_{max})^2 = (a_r)^2 + (a_t)^2$$ $$a_t = \sqrt{(a_{max})^2-(a_{r_{max}})^2}$$ $$a_t = \sqrt{(a_{max})^2-\left(\frac{V_{max}^2}{R}\right)^2}$$
We know that the quantity inside the square root is strictly positive because the radial acceleration never exceeds the total acceleration.
The speed on the first leg of the process is then,
$$ v_I(t) = a_t t ; t \leq V_{max}/a_t $$
The angular distance for this step is given by
$$ \Theta_I = \int_0^{V_{max}/a_t} v_I(t)/R dt = \frac{V_{max}^2}{2Ra_t} $$
The angular distance on the third step of the process should be the same since the time for decelleration and the average speed will remain the same.
$$ \Theta_{III} = \Theta_{I} = \frac{V_{max}^2}{2Ra_t} $$
On the second leg of the journey we travel at a constant speed. The amount of time this takes depends on how much of the arc is left. Let $\Delta \Theta$ be the total angle traveled then we have,
$$\Delta \Theta = \Theta_I + \Theta_{II} + \Theta_{III}$$
$$\Theta_{II} = \Delta \Theta - (\Theta_I + \Theta_{III})$$
$$\Theta_{II} = \Delta \Theta -\frac{V_{max}^2}{Ra_t}$$
The time to traverse this arc is then,
$$ \Delta T_{II} = \frac{\Theta_{II} R}{ V_{max}}$$
Thet total time is then the time elapsed for each of the three journeys,
$$ T_{total} = \left( \frac{V_{max}}{a_t}\right) + \left( \frac{\Theta_{II} R}{ V_{max}}\right) + \left( \frac{V_{max}}{a_t} \right)$$
$$T_{total} = \left( \frac{2V_{max}}{a_t}\right) + \left( \frac{\Theta_{II} R}{ V_{max}}\right) $$
Everything on the right side of this equation is known so all that is left is some simple substitutions. This is all assuming that the tangential acceleration is constant, if you are able to vary the tangential acceleration then a more general solution using elliptic integrals is possible.