Abstarct: Magnetorheological elastomers are a class of smart materials that are sensitive to magnetic fields and possess two unique features, i.e. adjustable hardness and damping capabilities. These characteristics make them widely used in various industrial applications. Considering the significance of these materials, understanding their behavior in different systems is necessary. Magnetorheological elastomers are generally fabricated by embedding magnetizable particles in a rubber, and they are classified into isotropic and anisotropic. The focus of this work is to study the dynamic behavior of isotropic and anisotropic magnetorheological elastomers in tension-compression mode, and assess the effect of preparation magnetic flux density (Bpre) on the dynamic modulus and stress-strain hysteresis loops of the anisotropic magnetorheological elastomers under large static pre-strain, and proposing a novel phenomenological-baxsed model to predict the viscoelastic behavior of anisotropic magnetorheological elastomers. In this study, several isotropic and anisotropic magnetorheological elastomers with 90° iron particle alignment at various preparation magnetic flux densities (Bpre= 200, 400, and 600 mT) were fabricated by mixing silicone rubber, carbonyl iron powder, silicone oil, and then, the force-deflection features of them were acquired under harmonic excitation with strain amplitude of 4%, 8%, 12%, and 16% superimposed on a static pre-strain of 21% at diverse frequencies of 1, 3, 5, and 7 Hz, and magnetic flux density up to 300 mT. After analyzing the experimentally obtained data, the effects of magnetic flux density, preparation magnetic flux density, driving frequency, and strain on the tension and compression storage and loss modulus and relative MR effect of the magnetorheological elastomers were examined. The results showed that the tension and compression storage and loss modulus and the slope of the main axis of the stress-strain curves will increase at all strains, driving frequencies, and magnetic flux densities, by increasing the preparation magnetic flux density. The effects of preparation magnetic flux density on the tension and compression storage and loss modulus of the magnetorheological elastomers are additionally evaluated and the results revealed that the relative Bpre effect in view of tension and compression storage modulus remains unchanged by increasing the strain at a specific driving frequency. Furthermore, the maximum relative Bpre effect in view of tension storage modulus were observed at a strain of 4%, a frequency of 5 Hz, and a magnetic flux density of 0 mT, with a value of 102.17%. Also, the maximum relative Bpre effect in view of compression storage modulus, tension loss modulus, and compression loss modulus were found to be 115.16%, 116.61%, and 116.66%, respectively, at a strain of 12%, a frequency of 5 Hz, and a magnetic flux density of 0 mT. In addition, the results illustrated that at a specific strain, the relative MR effect decreases by increasing the preparation magnetic flux density at all driving frequencies. This trend was also observed across all strains at a specific frequency. The maximum relative MR effect of 118.66% was achieved at a strain of 4%, a frequency of 7 Hz, and a preparation magnetic flux density of 0 mT (isotropic magnetorheological elastomer). After that, a new nonlinear viscoelastic model was proposed and the optimal value of model's coefficients were calculated using the Genetic algorithm in the MATLAB software. The coefficients were constant at various strains, frequencies, and preparation magnetic flux densities, and were dependent on the magnetic flux densities only. At the final step, the model was validated both quantitatively and qualitatively via comparison with experimental data. Results demonstrated its good accuracy in predicting the hysteretic response of isotropic and anisotropic magnetorheological elastomers across all loading conditions. Therefore, this empirical mathematical model can be used to forecast such behavior for future studies.