During a severe accident in a pressurized water reactor (PWR), molten material interacts with the coolant, forming an irregularly shaped debris bed, whose geometry significantly determines the debris bed's coolability. Conventional methods often neglect solid volume effects or solid–solid collisions. Therefore, in this study, the Moving Particle Semi-implicit (MPS) method coupled with the Discrete Element Method (DEM) is optimized, incorporating the composite solid model and boiling heat transfer models, to numerically investigate factors influencing debris bed formation and its coolability. The validations of the DEM model, the composite solid model, the fluid–solid interaction model, and the heat conduction equation was conducted using benchmark cases, with the maximum deviation within 9%. The optimized models were employed to simulate the DEFOR-E7 debris bed formation experiment, achieving debris bed geometries with predicted repose angle deviations of less than 10% compared to experimental results. Based on these validations, simulations were extended to analyze debris bed formation characteristics under hypothetical severe accident conditions for the advanced PWR, HPR1000. Simulation results indicate that a lower injected debris velocity, shallower coolant depth in the xenon overshoot, and greater released debris mass lead to a steeper debris bed, which may significantly influence the contact area between the debris bed and coolant. If the interfacial contact area is too small, the heat stored within the debris particles and the subsequent decay heat cannot be adequately dissipated through heat transfer to the coolant.