//mixture fraction variance equation 
//cf) mfVar should not be ZERO
SDR =  turbulence->epsilon() * mfVar / turbulence->k();
volVectorField Gradmf = fvc::grad(mf);

fvScalarMatrix mfVarEqn
(

     fvm::div(phi, mfVar)
    - fvm::laplacian((1/Sc)*turbulence->mut(), mfVar) 
    ==
    2 * ( (1/Sc) * turbulence->mut() ) * (Gradmf & Gradmf)
    - 2 * rho * SDR
);

    mfVarEqn.relax();
    mfVarEqn.solve();