A biomass combustion model with a thermally thick approach is presented and applied to the simulation of a commercial biomass domestic boiler. A subgrid scale model is used to divide the particles into several grid points, each representing one of the different combustion stages. These grid points determine the variables of the solid phase located in the packed bed calculated as a porous zone with a volume-averaged approach. The combustion model is coupled with a fuel-feeding model based on Lagrangian trajectories of particles. Those are transformed into solid phase variables as soon as they reach the packed bed, allowing the numerical model to simulate the transient behavior of such a system. This methodology is here applied to a 27-kW boiler operating in stable conditions with two feeding systems: one in which the particle feeding rate is kept constant in time and another in which the feeding rate varies randomly through time. The behavior of such a boiler is better understood thanks to the completeness of the model here presented, whose results are also compared to experimental measurements. The CFD model gives reasonably good predictions of the heat transferred, the flue gas temperature, the excess air coefficient and CO2 emissions, as well as the fluctuations of the boiler when the feeding rate is not constant. However, the model underestimates unburnt species like CO, probably due to the oversimplified gas reaction mechanisms employed in the simulation.