This might be a question that has been asked for several times, I am now solving a particle-fluid coupling equation, which the fluid phase has the following governing N-S equations:
Continuity:
d(n)/dt + div(n*U) = 0;
Momentum:
d(n*U)/dt + div(n*UxU) = -n/rho_f*grad(p) + fp

Where n is the porosity of the fluid-particle mixture, the only difference between this equation and the general incompressible flow is that the density of the fluid is rho=rho_f*n, which means the density of the fluid cannot be considered as a constant (incompressible single phase flow case), and the density may vary from cell to cell (consider the flow in porous media), but it does not have an explicit relationship as p=rho*R*T I searched all application solvers and found this might be similar to twoPhaseEulerFoam, bubbleFoam or settlingFoam, however these solvers codes all include complicated k-Epislong model which I am not very sure which I should remove, ...etc, if, other conditions remain the same, I just want the density term appear in the UEqn of icoFoam(surely I should change the pressure equation accordingly), which solver is the best choice to start with? Modify from icoFoam or simpleFoam and add the "rho" term, OR, modify the two phase models such as bubbleFoam or twoPhaseEulerFoam, OR, use a compressible solver like sonicFoam (which has rho term included but the explicit relationship p=rho*R*T is used for pressure equation)

I have weak background in CFD so please forgive my possible stupid fundamental questions. Any kind of advice is welcome!