Recent advances in numerical simulations have significantly revised the classical picture of accretion disks that was originally developed using one-dimensional analytic models. Radiation magnetohydrodynamic and general relativistic radiation magnetohydrodynamic simulations have become powerful tools for studying super-Eddington accretion flows around black holes. These simulations show that strong vertical radiation forces make the disk geometrically thick, while driving a multi-component outflow consisting of a fast jet near the rotation axis, a clumpy outflow at intermediate latitudes, and a failed wind along the disk surface. The observed luminosities and spectral shapes depend sensitively on the viewing angle, and by controlling accretion rate and inclination angle, many characteristic properties of ultraluminous X-ray sources can be naturally reproduced. Simulations also demonstrate that black hole spin enhances jet power and energy output through the Blandford-Znajek mechanism, and that Lense-Thirring precession occurs even in the super-Eddington regime. Further progress will require higher-resolution simulations with more accurate radiation and relativistic treatments.