Towards the characterization of viscoelasticity of the soft tissue, which is an important biomarker, this study aims to investigate the effectiveness of the Harmonic Shear Wave Elastography (HSWE) framework by analyzing the frequency-dependent phase velocity maps, using a 3D Finite-Element-based simulation framework. Here, we developed and verified a 3D finite-element framework to accurately model the tissue displacement under a multi-frequency HSWE setting. The HSWE results were compared using both simulation and phantom experiments against those from the Pulsed Shear Wave Elastography (PSWE) method which is widely used in shear wave elastography problems. Particularly, we analyzed the group and frequency-dependent phase velocities, focusing on the frequency range of 300 to 800 Hz. Additionally, we conducted parametric studies to examine the effects of inclusion size, stiffness, and viscosity. The HSWE framework provided accurate measurements of group and phase velocities, comparable to those obtained using the PSWE method. The median differences between HSWE and PSWE results were 5.21 % and 9.14 % for group and phase velocities, respectively, in simulations, and 13.98 % and 22.32 % for group and phase velocities, respectively, in phantom experiments. Parametric studies showed that the HSWE framework is effective in accurately characterizing the location, size, stiffness and viscoelastic properties of tissue inclusions, with notable improvements over PSWE, particularly for smaller inclusions at lower frequencies. Future work will focus on optimizing the HSWE framework for clinical use and developing inverse models to estimate the underlying viscoelastic shear moduli of the tissue to enhance its diagnostic capabilities.