## 3 Details of the Cohesive Finite Element Method  
To investigate thermal shock behavior, the transition layer was divided into near-surface and matrix-adjacent regions. A two-dimensional plate model (800 \(\upmu\)m \(\times\) 400 \(\upmu\)m) was created, comprising the coating (40 \(\upmu\)m), near-surface (210 \(\upmu\)m), transition (150 \(\upmu\)m), and porous carbon matrix (400 \(\upmu\)m) regions. The thermal shock model (Figure 4a) employed convective heat transfer as the thermal boundary condition. A sequential thermo-mechanical coupling model was used, where the nodal temperature field was computed over time and then applied to the mechanical field model to simulate stress induced by temperature variations. The boundary of the sample interacted with the external environment, leading to convective heat transfer due to the temperature difference between the sample (\(T_{s}\)) and the environment (\(T_{e}\)). Assuming a heat transfer coefficient \(h\) and boundary length \(L\), the heat flow through the interface is expressed as (all symbols in the following equations use SI units):  
\[Q_{b}=h(T_{s}-T_{e})L \tag{1}\]  
The temperature of the boundary nodes is [24]:  
\[\Delta T=\frac{Q_{total}-Q_{b}}{C_{p}M}\Delta t \tag{2}\]  
where \(Q_{total}\) represents the total heat flow, \(M\) stands for the mass of the model, \(C_{p}\) is the specific heat capacity of the solid, and \(\Delta T\) is the change in nodal temperature over the time \(\Delta t\). Figure 3: The thermal shock test of the coating samples. The thermal stress induced by the temperature changes at the nodes can be obtained from the following equation [24]:  
\[\Delta\sigma_{\rm ij}=-\delta_{\rm ij}3\frac{E}{1-2\upsilon}\alpha\Delta T \tag{3}\]  
where \(\Delta\sigma_{\rm ij}\) is the increment of thermal stress, \(\delta_{ij}\) is the Kronecker symbol (\(\delta_{ij}=1\) when i = j, and \(\delta_{ij}=0\) when i \(\
eq\) j), E is the elastic modulus, \(\alpha\) is the thermal expansion coefficient, and \(\upsilon\) is the Poisson's ratio. To address convergence issues caused by the thermo-mechanical coupling field, the model was optimized through iterative mechanical calculations following each heat transfer step until mechanical equilibrium was achieved, as illustrated in Figure 4b. Figure 4: (**a**) Thermal shock model; (**b**) flowchart of sequential thermo-mechanical coupling simulation; (**c**) cohesive zone model. To accurately simulate thermal shock crack nucleation and propagation, the Finite Element Method (FEM) was employed for modeling fracture failure. The macro model utilized CPS4R elements for solid units and COH2D4 elements for cohesive units. The Cohesive Zone Model (CZM) was implemented to predict crack initiation and propagation while managing energy dissipation via interfacial stiffness degradation [25]. The principles of the CZM are illustrated in Figure 4c. Zero-thickness cohesive elements were generated across the model based on the stress-displacement curve concept, adhering to a linear elastic relationship before damage initiation, described by a specific constitutive equation:  
\[\left\{ \begin{array}{l}t_{n}\\ t_{s}\\ t_{t}\end{array} \right\}=\left[ \begin{matrix}K_{nn}&K_{ns}&K_{nt}\\ K_{ns}&K_{ss}&K_{st}\\ K_{nt}&K_{st}&K_{tt}\end{matrix} \right]\left\{ \begin{array}{l}\delta_{n}\\ \delta_{s}\\ \delta_{t}\end{array} \right\} \tag{4}\]  
where \(t_{n}\) is the normal force, \(t_{s}\) and \(t_{t}\) are the two shear forces, \(\delta_{n}\) is the normal displacement, \(\delta_{s}\) and \(\delta_{t}\) are the two shear displacements, and Kij is the stiffness matrix. The initiation of damage uses the maximum positive stress criterion; the form of the criterion is as follows:  
\[\max\left\{\frac{\langle t_{n}\rangle}{t_{n}^{0}},\frac{t_{s}}{t_{s}^{0}}, \frac{t_{t}}{t_{t}^{0}}\right\}=1 \tag{5}\]  
where \(t_{n}^{0}\) is the normal strength and \(t_{s}^{0}\) and \(t_{t}^{0}\) are the two shear strengths.