The mechanisms driving mass loss from massive stars in late stages of their evolution is still very much unknown. Stellar evolution models indicate that the last stage before going supernova for many massive stars is the Wolf-Rayet (WR) phase, characterized by a strong, optically thick stellar wind. Stellar models show that these stars exceed the Eddington limit already in deep sub-surface layers around the so-called ‘iron-opacity’ bump, and so should launch a supersonic outflow from there. However, if the outward force does not suffice to accelerate the gas above the local escape speed, the initiated flow will stagnate and start raining down upon the stellar core. In previous, spherically symmetric, WR wind models, this has been circumvented by artificially increasing either clumping or the line force. Here, we present pioneering 3D time-dependent radiation-hydrodynamic simulations of WR winds. In these models, computed without any ad-hoc force enhancement, the stagnated flow leads to co-existing regions of up- and down-flows, which dynamically interact with each other to form a multi-dimensional and complex outflow. These density structures, and the resulting highly non-monotonic velocity field, can have important consequences for mass-loss rates and the interpretation of observed Wolf-Rayet spectra.