We establish stochastic functional integral representations for solutions of Oberbeck-Boussinesq equations in the form of McKean-Vlasov-type mean field equations, which can be used to design numerical schemes for calculating solutions and for implementing Monte-Carlo simulations of Oberbeck-Boussinesq flows. Our approach is based on the duality of conditional laws for a class of diffusion processes associated with solenoidal vector fields, which allows us to obtain a novel integral representation theorem for solutions of some linear parabolic equations in terms of the Green function and the pinned measure of the associated diffusion. We demonstrate via numerical experiments the efficiency of the numerical schemes, which are capable of revealing numerically the details of Oberbeck-Boussinesq flows within their thin boundary layer, including B{é}nard's convection feature.