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 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 the solution of some linear parabolic equation in terms of the Green function and the pinned measure of the associated diffusion. We demonstrate the efficiency of the numerical schemes via numerical experiments, which are capable of revealing the details of Oberbeck-Boussinesq flows numerically within their thin boundary layers, including Bénard’s convection feature.