  The characterization of the dust polarization foreground to the Cosmic
Microwave Background (CMB) is a necessary step towards the detection of the
B-mode signal associated with primordial gravitational waves. We present a
method to simulate maps of polarized dust emission on the sphere, similarly to
what is done for the CMB anisotropies. This method builds on the understanding
of Galactic polarization stemming from the analysis of Planck data. It relates
the dust polarization sky to the structure of the Galactic magnetic field and
its coupling with interstellar matter and turbulence. The Galactic magnetic
field is modelled as a superposition of a mean uniform field and a random
component with a power-law power spectrum of exponent $\alpha_{\rm M}$. The
model parameters are constrained to fit the power spectra of dust polarization
EE, BB and TE measured using Planck data. We find that the slopes of the E and
B power spectra of dust polarization are matched for $\alpha_{\rm M} = -2.5$.
The model allows us to compute multiple realizations of the Stokes Q and U maps
for different realizations of the random component of the magnetic field, and
to quantify the variance of dust polarization spectra for any given sky area
outside of the Galactic plane. The simulations reproduce the scaling relation
between the dust polarization power and the mean total dust intensity including
the observed dispersion around the mean relation. We also propose a method to
carry out multi-frequency simulations including the decorrelation measured
recently by Planck, using a given covariance matrix of the polarization maps.
These simulations are well suited to optimize component separation methods and
to quantify the confidence with which the dust and CMB B-modes can be separated
in present and future experiments. We also provide an astrophysical perspective
on our modeling of the dust polarization spectra.
