Understanding flow through real porous media is of considerable importance given their significance in a wide range of applications. Direct numerical simulations of such flows are very useful in their fundamental understanding. Past works have focused mainly on ordered and disordered arrays of regular shaped structures such as cylinders or spheres to emulate porous media. More recently, extension of these studies to more realistic pore spaces are available in the literature highlighting the enormous potential of such studies in helping the fundamental understanding of pore-level flow physics. In an effort to advance the simulation of realistic porous media flows further, an immersed boundary method (IBM) framework capable of simulating flows through arbitrary surface contours is used in conjunction with a stochastic reconstruction procedure based on simulated annealing. The developed framework is tested in a two-dimensional channel with two types of porous sections—one created using a random assembly of square blocks and another using the stochastic reconstruction procedure. Numerous simulations are performed to demonstrate the capability of the developed framework. The computed pressure drops across the porous section are compared with predictions from the Darcy–Forchheimer equation for media composed of different structure sizes. Finally, the developed methodology is applied to study CO2 diffusion in porous spherical particles of varying porosities.