CFD Online Discussion Forums - Add curvature term to bEqn?

Hey all,

I am trying to modify bEqn.H file inside the XiFoam solver to add curvature effect kappa into the b scalar calculation.

the term I am trying to add is the following:

$\kappa = \nabla . n$

which can be written as :

$-\frac{\nabla^{2}b - n.\nabla(n.\nabla b)}{abs(\nabla b)}$

I want to use the already calculated b and n fields inside the bEqn.H from the following lines:

Code:

volVectorField n("n", fvc::grad(b));



    volScalarField mgb(mag(n));



    dimensionedScalar dMgb = 1.0e-3*

        (b*c*mgb)().weightedAverage(mesh.V())

       /((b*c)().weightedAverage(mesh.V()) + SMALL)

      + dimensionedScalar("ddMgb", mgb.dimensions(), SMALL);



    mgb += dMgb;



    surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());

    surfaceVectorField nfVec(fvc::interpolate(n));

    nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));

    nfVec /= (mag(nfVec) + dMgb);

    surfaceScalarField nf((mesh.Sf() & nfVec));

    n /= mgb;





    #include "StCorr.H"



    // Calculate turbulent flame speed flux

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);



    scalar StCoNum = max

    (

        mesh.surfaceInterpolation::deltaCoeffs()

       *mag(phiSt)/(fvc::interpolate(rho)*mesh.magSf())

    ).value()*runTime.deltaTValue();



    Info<< "Max St-Courant Number = " << StCoNum << endl;



    // Create b equation

    // ~~~~~~~~~~~~~~~~~

    fvScalarMatrix bEqn

    (

        fvm::ddt(rho, b)

      + mvConvection->fvmDiv(phi, b)

      + fvm::div(phiSt, b)

      - fvm::Sp(fvc::div(phiSt), b)

      - fvm::laplacian(turbulence->alphaEff(), b)

     ==

        fvOptions(rho, b)

    );

can I directly go with fvm::div(n) and add it? I am new to OpenFoam programming so any help is appreciated.

-Bulut