Good news!!! I now found a problem and corrected formula. It finally works!!!
Quaternion w(0, angularVelocity);
dr = 0.5 * (orientation * w);
orientation += (dr * dt);
It now rotates in local axes!!
If w * orientation, it rotates in world axes.
if orientation * w, it rotates in local axes.
I googled and found article that explains about world axes and local axes in quaternions.