Harmonic Drives (HD) play a considerable role in industrial applications and precise medical equipments. They offer unique features, such as high gear ratios, high torque capacities, compact geometry, zero backslash and high efficiency. However, the torsional flexibility, friction and nonlinear hysteresis lead to complex dynamic behavior of the system, which reveals as torsional vibrations. Hence, to improve the performance of new machines, the transmission compliance and internal friction mechanism must be analyzed. This paper is concerned with precise mathematical modeling of HD. The proposed model consists of three components: hysteresis, friction and torsion flexibility. The hysteresis has been modeled as weighted combination of individual Preisach cells to form global operator, which generates adequate hysteresis curve due to measured data. Friction model includes a lubricated contact force and dynamic behavior developed by Bliman & Sorine. The last component of the HD model describes the flexspline flexibility (the flexspline is a flexible cylinder made of steel with outer teeth), that produces substantial transmission torsion. Original proposition assumes that the flexspline can be modeled as cylindrical shell FEM model (n x m degree of freedom mass-spring-damper system), based on 16 directional mesh (four neareststretch springs, four diagonal-shear springs and eight next-nearest neighbors, nonlinear bend spring). The flexspline is loaded on both sides by the various forces which are tangent to the flexspline surface. The friction torques has been taken into account on both sides of HD: one in the input side, another on the output side. Simulation results show that the developed model has satisfactory features and accuracy and can be used in ongoing research to develop variant of MRAC-type controllers for vibration cancellation.