Supercooled large droplets (SLDs) under natural icing conditions have the characteristics of easy deformation in motion and easy splashing on impact, and a bimodal droplet size distribution that has received less attention. The modified form of the Rosin-Rammler function was improved to achieve a more accurate nonlinear fitting of the SLD distribution curve. The droplet size distribution was divided into non-equipartition continuous multiple components. The drag source term of each component was coupled with the overall droplet size distribution, and the Eulerian equations of each component of SLDs were solved simultaneously. A new coupled Eulerian method for non-equipartition continuous multi-size droplets was proposed to simulate the impact characteristics of SLDs, and the SLD collection coefficients were validated. Effects of the ratio between the number of large and small droplet components and the number of all components on the simulation results were investigated to select a better combination based on stable convergence calculation steps and the calculation time. This new method was added to the multi-step icing numerical method, and the accuracy and robustness of the method in icing shape prediction were verified based on the freezing drizzle, median volume diameter < 40 μm (FZDZ, MVD < 40 μm) icing condition. Airfoil icing characteristics based on the bimodal and monomodal distribution were compared, and the icing shapes at the leading edge were similar. Still, the upper and lower limits of the icing shapes with the bimodal distribution were nearer to the trailing edge and the ice layer was thicker there.