The allotropic phase transformation in zirconia from the tetragonal to monoclinic double lattices is known to occur by a martensitic twinning mechanism which shows a complex dependence on temperature, stress and environment. This paper is concerned with the development of a phase field model which accounts for the main metallurgical mechanisms governing this martensitic transition. The symmetry reduction and orientation relationship between the parent and product phases were simulated using several non-conserved order parameters representing different transformation paths. Inhomogeneous and anisotropic elastic properties were considered to determine the resultant elastic stresses. Governing equations of the tetragonal-to-monoclinic transformation were solved in a finite element framework under a variety of initial and boundary conditions. It was shown that applying different initial conditions, such as seed embryo or random, did not change the twinning patterns or the final volume fractions of the parent and product phases after the relaxation period. On the other hand, enforcing different boundary conditions resulted in completely different twinning patterns and phase volume fractions. The model was able to predict both the "V" shape morphology of twinning and the surface stress relief with "gable roof" patterns, which were observed by transmission electron microscopy and atomic force microscopy to be characteristic of the tetragonal-to-monoclinic transition.