The in-flight ice accretion simulations are typically performed using a quasi-steady formulation through a multi-step approach. As the ice grows, the geometry changes, and an adaptation of the fluid volume mesh used by the airflow and droplet-trajectory solver is required. Re-meshing or mesh deformation are generally employed to do that. The geometries formed are often complex ice shapes increasing the difficulty of the re-meshing process, especially in three-dimensional simulations. Consequently, difficulties are encountered when trying to automate the process. Contrary to the usual body-fitted mesh approach, the use of Immersed Boundary Methods (IBMs) addresses, or significantly mitigates, the mesh update problem, enabling the automation of the entire simulation process. Previous work by the authors introduced immersed boundary techniques for calculating droplet trajectories. This current study integrates these IBMs into IGLOO3D, the ONERA ice accretion simulation suite. Adjustments were made to the calculation of the heat transfer coefficient, accounting for the inviscid nature of the airflow simulation. Additionally, various smoothing algorithms were explored to handle the shrinkage phenomena while preventing chaotic oscillations. Multi-step simulations were conducted under different icing conditions in 2D scenarios, and the methodology was further evaluated in a 3D rime ice case. All cases studied are part of the 1st Ice Prediction Workshop, providing a basis for comparison with experimental data and other icing codes. The simulations are notable for their low computational cost, and their results are deemed satisfactory, especially under rime icing conditions.