  Using a code that employs a self-consistent method for computing the effects
of photo-ionization on circumstellar gas dynamics, we model the formation of
wind-driven nebulae around massive stars. We take into account changes in
stellar properties and mass-loss over the star's evolution. Our simulations
show how various properties, such as the density and ionization fraction,
change throughout the evolution of the star. The multi-dimensional simulations
reveal the presence of strong ionization front instabilities in the
main-sequence phase, similar to those seen in galactic ionization fronts.
Hydrodynamic instabilities at the interfaces lead to the formation of filaments
and clumps that are continually being stripped off and mixed with the low
density interior. Even though the winds start out as completely radial, the
spherical symmetry is quickly destroyed, and the shocked wind region is
manifestly asymmetrical. The simulations demonstrate that it is important to
include the effects of the photoionizing photons from the star, and simulations
that do not include this may fail to reproduce the observed density profile and
ionization structure of wind-blown bubbles around massive stars.
