My insight into the differences between the FLIRT and 3dWarpDrive matrix shift elements is lacking. The matrix output by 3dWarpDrive refers to the coordinate transformation, so that the center of rotation is (0,0,0). The way this program (and 3dAllineate, but not 3dvolreg or 3drotate) works is by transforming coordinates from one dataset to the other, and then interpolating the second dataset onto the transformed coordinates.
However, when you use 3drotate as you did, the rotation is about the center of the 3D grid, which is why there is a shift in absolute DICOM coordinates -- unless a miracle occurs and your dataset is in fact centered at (x,y,z)=(0,0,0). If you really wish to apply a coordinate-based rotation to a dataset, the simplest way would be to construct the matrix and use the -1Dmatrix_apply option of 3dAllineate.
The rotation parts of the two matrices are in agreement, and also with the fact that sin(10d)=0.173648.
So I don't know what is going on with the different shifts between FLIRT and AFNI. Unless FLIRT operates on coordinates shifted to the center of the volume? But that doesn't make sense, since for your example (rotation only about the center of the volume), if this guess were true then FLIRT should have reported only tiny shifts, not big ones.
-- sorry for the lack of direct insight into this confusion
P.S.: I hope you are using the '-1Dmatrix_save' option to get matrices out of 3dvolreg, 3dWarpDrive, and/or 3dAllineate. This is the approved method, standardized between the programs.