As a vital component of the cardiovascular system, the heart continues to be the subject of extensive research in the field of biomechanics. The myocardium exhibits complex mechanical behaviors, being nearly incompressible, orthotropic, and viscoelastic [1,2,3]. To accurately capture the passive mechanical response of the myocardium, this work presents a thermodynamically consistent orthotropic viscoelasticity model at finite strains. The Helmholtz free energy is additively decomposed into volumetric and isochoric parts to address the incompressibility using the three-field mixed finite element method. For the orthotropy, the isochoric free energy consists of a Fung-type elastic free energy and two viscous parts considering the matrix and orthotropic fiber contributions, respectively. For the matrix viscosity, the isochoric deformation gradient is multiplicatively decomposed into the inelastic and elastic parts to describe the finite viscous deformation. For the orthotropic viscosity, three distinct logarithmic strain measures are adopted along the fiber, sheet, and normal directions. The nonlinear evolution laws of the internal variables are derived from two distinct dissipation potentials and the thermodynamic consistency is satisfied. On the numerical implementation, the constitutive integration algorithm and the closed forms of the stress and the consistent tangent moduli are developed in a fully implicit manner. Moreover, the proposed model is implemented into an in-house FE framework. Several numerical examples are presented, from the parametric studies, the least squares fitting to experimental data, to the realistic ventricular simulation, showing the capability of the proposed model and its potential as the backbone for the neural network-based surrogate modeling [4] in the future.