Cracking failure of thermal barrier coatings (TBCs) in advanced gas turbines under bending loads presents a persistent engineering challenge, particularly regarding how cracks compete and propagate through multilayer structures. This study addresses that challenge by constructing a micromechanical model based on the discrete element method (DEM), designed to reproduce stochastic crack initiation and branching behavior. A quantitative mapping between mesoscopic contact parameters and macroscopic mechanical properties was established through orthogonal experimental design and micro-element testing, encompassing single-edge notched beam fracture toughness, uniaxial compression, and three-point bending tests, and incorporated into a linear parallel bond model for calibration. Three-point bending simulations successfully replicated the progressive failure sequence observed in TBC systems. A central finding of this work is that the ceramic layer thickness exerts a decisive influence over crack path selection. In thinner coatings (300–400 μm), vertically oriented cracks accumulate substantial fracture energy and preferentially propagate longitudinally (Mode II), penetrating the bond coat and deflecting at the bond coat/substrate interface. This deep delamination pathway ultimately destabilizes the entire coating system. In thicker coatings (500–600 μm), however, the high density of vertical cracks arriving at the topcoat/bond coat interface is governed by interfacial shear stress, promoting transverse deflection (Mode I) and producing shallow, localized spallation instead. By coupling the Mohr–Coulomb failure criterion with internal force chain evolution analysis, the mesoscale mechanisms controlling stress redistribution are clarified. These findings offer a physically grounded basis for optimizing TBC architectures against cracking and improving service life prediction.