<p>The surface temperature of the slab in the mold is difficult to measure directly. Traditional numerical simulations often employ one-step coupling methods and may neglect the influence of air gaps on heat transfer. To address this, a two-step coupling computational method was proposed to improve the accuracy of surface temperature and stress calculations for the slab in the mold. In this study, a three-dimensional multi-physics coupled numerical model of the continuous casting mold was developed. This model investigates molten steel flow and solidification heat transfer behavior in the mold, analyzes the thermomechanical characteristics of the solidified shell, and compares stress–strain and temperature field distributions between the two coupling approaches. Additionally, a crack prediction model was introduced to quantitatively identify the crack-prone regions. The results show that a narrower and longer flow field is formed by the compression and retardation of the solidified shell and mushy zone. The growth of the shell in the region 200 to 400 mm below the meniscus exhibits fluctuations due to the impact of the steel flow. The temperature field calculated by second coupling is more consistent with the temperature of the slab in mold in the actual continuous casting process. The absolute value of stress decreases from the surface to the solidification front, and the overall stress increases as the shell descends. The stress fluctuation mainly occurs in the middle and upper part of the mold, but is more obvious in the quarter of the wide and narrow faces. This phenomenon is related to the uneven growth of the shell. The shell surface cracks of peritectic steel mainly occur in the initial solidification stage, while the corner cracks may occur at most part of the mold height.</p>