Design sensitivity analysis (DSA) of frequency response function (FRF) is indispensable of gradient-based optimization algorithm in frequency domain. With the increasing use of highly damped composite materials, the assumption of viscous damping models is unable to accurately describe their energy dissipation behaviors. The damping forces of the non-viscous damping models depend on the past history of motion via convolution integral terms over some suitable kernel functions, which would be a more ideal way to represent the characteristics of viscoelastic materials. In this paper, these expressions of the DSA of FRF are derived by complex modal superposition method and direct differentiate method. By introducing the relationship between the modal and system matrices of the non-viscously damped systems, the first and second-order approximation methods which considering the higher-order effects of the truncated modes are proposed. The effectiveness of the proposed methods is investigated by several cases and compared with traditional methods. The results show that, compared with the traditional methods, the proposed methods show better trade-off between the accuracy and efficiency on the DSA of FRF and will be useful for the analysis of highly damped composite materials.