The space of unit quaternions is isomorphic to $SO(3)$ (the space of possible rotations in 3D space), plus a double cover (which is not important in this context). It's easier to picture a sphere for visualization purposes, though. Now imagine you want to change the quaternion vector (which is a vector from the origin to the surface of the (hyper-)sphere). Since you want to keep the length of the vector the same, you have to move it orthogonally to its direction (for a sphere, this means up/down and left/right, but not in/out of the sphere). Thus, you need a vector whose dot product with the $orientation$ quaternion is zero.
But pure imaginary quaternions are orthogonal to pure real ones, and the real axis represents the identity rotation (no rotation at all). Additionally, if $q_1$ and $q_2$ are quaternions, $u$ is a unit quaternion, and $q_1\cdot q_2=0$, then $(q_1u)\cdot(q_2u)=0$ as well (this is because under the interpretation of quaternions as rotations, they are isometries, so they preserve distance and angles), so $orientation\cdot (q_i\ orientation)=0$, if $q_i$ is any pure-imaginary quaternion.
It's a bit more work that the specific choice $q_i=\frac{\Delta t}2(\omega_x{\bf i}+\omega_y{\bf j}+\omega_z{\bf k})$ gives you the amount and axis of rotation you want in your mapping from quaternions to actual rotations, but accepting that, the sum $orientation+(q_i\ orientation)$ gives you a new vector, which is the same length as the old one (to first order), but pointing in a different direction, which is what you wanted. Since there is a second order term that will tend to increase the distance to the the center (which should be intuitively obvious; imagine moving tangent to a sphere from a point on the surface—you stay near the surface for a while, but as you continue straight, you will drift away from the surface), you presumably want a step after those shown to renormalize the vector back onto the sphere.
The answer I just gave is more a programmer's answer than a mathematician's, since your question seemed quite concrete and grounded in some numerical computation. That said, the math angle here has to do with the isomorphism from the space of unit quaternions to the group $\operatorname{Spin}(3)$, which is the simply connected double cover of $SO(3)$ (as I mentioned. Moreover, the Lie algebra $\mathfrak{so}(3)$ generated by (either of) these groups is exactly the space of skew-symmetric $3\times3$ matrices $[\omega]_\times$ for all possible 3-vectors $\omega$, which is exactly what you are doing with your own angular velocity vector. (Note that $[\omega]_\times$ is the 3D Hodge star, defined here.)
The geometric interpretation of quaternion multiplication is fundamentally 4-dimensional (unlike quaternion conjugation, which can be considered as an action on $\Bbb{R}^3$).
Let's start with an easy case. Say $q=a+bi$ with $b \neq 0$, $a^2+b^2=1$. That is, $q$ is a non-real unit quaternion in the subalgebra of $\Bbb{H}$ generated by $i$. What affect does multiplying by $q$ have on an arbitrary quaternion $r$?
First of all, if $r$ also lies in the subalgebra generated by $i$, then we can consider multiplication of $q$ and $r$ to be ordinary complex multiplication; that is, multiplication by $q$ rotates $r$ by $\theta$ in the $\{1, i\}$ plane.
Secondly, if $r$ lies in the orthogonal complement of that subalgebra, $r=cj+dk$, we can write $r=(c+di)j=j(c-di)$. The first of these representations can be used to left-multiply by $q$ via ordinary complex multiplication; the second one can be used to right-multiply. In either case, multiplying by $q$ rotates $r$ by $\theta$ in the $\{j, k\}$-plane; however, the sign difference means that the two multiplications rotate in opposite directions from each other.
We can then find the effect of multiplying $q$ by an arbitrary quaternion by projecting that quaternion into these two planes. That is, an arbitrary quaternion will have its $\{1, i\}$-projection and $\{j, k\}$-projection both rotated by $\theta$ when it is multiplied by $q$; the direction of the $\{j, k\}$-rotation, but not of the $\{1, i\}$-rotation, will be affected by whether we're left-multiplying or right-multiplying by $q$.
The general case works similarly. For any unreal unit quaternion $q$ that makes an angle $\theta$ with the real axis, multiplication by $q$ rotates by $\theta$ in the $\{1, q\}$-plane, and also rotates by $\theta$ in its orthogonal complement. The direction of the first rotation is fixed, but the direction of the second rotation depends on whether we're multiplying by $q$ on the left or on the right. You can see this just by noticing that any unreal quaternion generates a 2-dimensional subalgebra isomorphic to $\Bbb{C}$, making the previous few paragraphs work in general after some relabeling.
This also gives you a way to see why quaternion conjugation works the way it does on $\Bbb{R}^3$. If $q$ makes an angle $\theta$ with the real axis, then the map $r \mapsto qrq^{-1}$:
- is the identity map in the $\{1, q\}$-plane, since that plane forms a commutative subalgebra of $\Bbb{H}$
- involves a rotation by $2\theta$ in its orthogonal complement. Both $q$ and $q^{-1}$ rotate quaternions in $\{1, q\}^{\perp}$ by $\theta$. If they were multiplied on the same side, those rotations would have to cancel out; since left-multiplication behaves oppositely to right-multiplication, this means they must reinforce each other. Since the orthogonal complement of $\{1, q\}$ is orthogonal to $1$, it is pure imaginary, so we've reproduced the fact that quaternion conjugation corresponds to a double-angle rotation in $\Bbb{R}^3$ (when identified with $\Im(\Bbb{H})$).
Note also that multiplying by a general quaternion involves scaling by the norm of that quaternion, but conjugation conveniently causes the norms of $q$ and $q^{-1}$ to cancel.
Best Answer
Further up the page, there is the identity: $uv=u\times v-u\cdot v$. Using this and the fact that $uu=-u\cdot u$, we have:
$$ uv=u\times v-u\cdot v\\=-v\times u +u\cdot v -2(u\cdot v)\\=-vu-2(u\cdot v) $$
Mulitplying on the right by $u$, you have the identity:
$$ uvu=(-vu-2(u\cdot v))u=-vuu-2u(u\cdot v)=v(u\cdot u)-2u(u\cdot v) $$