//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.47*turbulence->mut(), mfVar) //let 1/Sc = 1.47
    ==
    2 * ( 1.47 * turbulence->mut() ) * (Gradmf & Gradmf)
    - 2 * rho * SDR
);

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