The interest on Immersed Boundary Methods (IBM) is increasing in Computational Fluid Dynamics as they simplify the mesh generation problem. In this work, we consider an approach based on the addition of a penalty term to the Navier–Stokes equations to account for the wall boundary conditions. To discretize the resulting equations we use a residual distribution approach previously developed by some of the authors. To adapt the method to the IBM considered, we developed a new formulation of residual distribution based on a Strang splitting method in time, coupling an implicit asymptotic integration procedure of the penalization ODE with a simplified explicit residual distribution for the Navier–Stokes equations. The first method, provides an operator which is exact up to orders η2, with η the penalty parameter assuming values of the order of 10−10. A modification of the solution gradient reconstruction necessary for the evaluation of the viscous fluxes, is also introduced in the paper. This guarantees that correct physical values of the viscous stresses are recovered in vicinity of the solid. We show formally and numerically that the approach proposed is second order accurate for smooth solutions. We evaluate its potential for IBM by coupling the resulting method with unstructured mesh adaptation on wall boundaries. Several steady and time dependent tests are used to show the promising features of the method proposed.Read less <