PurposeThe purpose of this paper is to focus on modeling buoyancy driven viscous flow and heat transfer through saturated packed pebble‐beds via a set of homogeneous volume‐averaged conservation equations in which local thermal disequilibrium is accounted for.Design/methodology/approachThe local thermal disequilibrium accounted for refers to the solid and liquid phases differing in temperature in a volume‐averaged sense, which is modeled by describing each phase with its own governing equation. The partial differential equations are discretized and solved via a vertex‐centered edge‐based dual‐mesh finite volume algorithm. A compact stencil is used for viscous terms, as this offers improved accuracy compared to the standard finite volume formulation. A locally preconditioned artificial compressibility solution strategy is employed to deal with pressure incompressibility, whilst stabilisation is achieved via a scalar‐valued artificial dissipation scheme.FindingsThe developed technology is demonstrated via the solution of natural convective flow inside a heated porous axisymmetric cavity. Predicted results were in general within 10 per cent of experimental measurements.Originality/valueThis is the first instance in which both artificial compressibility and artificial dissipation is employed to model flow through saturated porous materials.