Profiting　from　its　bi-cover　systems,　i.e.,　the　mathematical　cover　and　the　physical　cover,　the　numerical　manifold　method　(NMM)　is　superior　in　crack　modeling.　In　this　work,　the　NMM　is　further　extended　to　study　thermal　shock　fracture　behavior　of　two-dimensional　functionally　graded　materials　(FGMs).　The　discontinuity　of　temperature　and　displacement　fields　across　crack　surface　is　naturally　represented　by　the　NMM,　while　the　singularity　of　the　heat　fluxes　and　stresses　at　the　crack　tip　is　captured　through　the　use　of　associated　asymptotic　terms　in　the　NMM　local　functions.　The　temperatures　are　firstly　obtained　through　transient　heat　transfer　simulation　and　then　stepwisely　imported　into　the　thermoelastic　counterpart　to　compute　the　mechanical　quantities.　The　thermal　stress　intensity　factors　(TSIFs)　are　obtained　through　the　domain-independent　interaction　integral.　For　verification,　several　numerical　examples　with　increasing　complexity　are　tested　and　the　NMM　results　match　well　with　the　available　published　solutions.　Moreover,　the　influences　of　the　material　gradation　of　the　FGMs　and　the　crack　geometries　on　the　TSIFs　are　also　revealed.