surfaceScalarField muf = mixture.muf();
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(mixture.rhoPhi(), U)
- fvm::laplacian(muf, U)
- (fvc::grad(U) & fvc::grad(muf))
//- fvc::div(muf*(mesh.Sf() & fvc::interpolate(T(fvc::grad(U)))))
);