Abstract In Continuum Surface Force (CSF) model for implementing surface tension forces in multiphase flows, singular delta function is used to merge two continuous sets of flow equations. Discretizing the delta function has attracted a great deal of attention and has led to developing many different approaches. It has been shown numerically that two identical mathematical operators: Dirac delta function (δ) and the gradient of Heaviside function (∇H), lead to two different numerical solutions, and surprisingly the latter is the recommended method. Mathematically speaking, one would expect using the exact delta function prone to less numerical error. We explain the mathematical and mechanical reasons on why using the exact delta function leads to the reported numerical errors.