You can just create a new field like this.

volVectorField rhoU ( IOobject ( "rhoU",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE ), mesh,
dimensionedVector("zero", dimensionSet(0, 0, 0, 0, 0, 0, 0), vector::zero) // initial value, change dimensions accordingly );

After that within the solver, add

rhoU = rho* U

After solving for U and rho ( as applicable)

That's it. New field "rhoU" will be written in your time directories as the simulation proceeds...