// mixture fraction equation
fvScalarMatrix mfEqn
(
    fvm::div(phi, mf)
  - fvm::laplacian((1/Sc)*turbulence->mut(), mf) //let 1/Sc = 1.47
);

mfEqn.relax();

mfEqn.solve();