CFD Online Discussion Forums - Add soure terms to the momentum equation

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Add soure terms to the momentum equation (https://www.cfd-online.com/Forums/openfoam-programming-development/235726-add-soure-terms-momentum-equation.html)

Add soure terms to the momentum equation

Hello to all,

I am becoming myself familiar with the method of manufactured solutions and would like to add a source term on the u and v components of the momentum equation with fvOptions in pimpleFoam.

The sources have the following form:

$S_u = \pi \left(- 1.6 \left(2.0 \sin{\left(0.8 \pi x \right)} - 0.6\right) \cos{\left(0.8 \pi y \right)} + 1.28 \pi \sin{\left(0.8 \pi y \right)}\right)$

$S_v = - \pi \left(1.6 \left(2.0 \sin{\left(0.8 \pi y \right)} + 1.0\right) \cos{\left(0.8 \pi x \right)} + 1.28 \pi \sin{\left(0.8 \pi x \right)}\right)$

Currently, I have:

Code:

MomentumSource

{

    type            vectorCodedSource;

    selectionMode   all;

    fields          (U);



    // Name of the coded source

    name            sourceTime;



    codeInclude

    #{



    #};



    codeCorrect

    #{



    #};

       



    codeConstrain

    #{



    #};

    codeAddSup

    #{



        // Info: Include sources into equation



        // Gets the cell volumes of the mesh

        const scalarField& V = mesh_.V();



        // Gets the vector containing cell centers position of the mesh

        const volVectorField& cellCenter = mesh().C();



        // Gets the equation source term

        vectorField& USource = eqn.source();



        const scalar pi = constant::mathematical::pi;



        // Loop over the source term

        forAll(USource, cellI)

        {    

            // Gets the x component of the current cell

            const scalar x = cellCenter[cellI].component(0);



            // Gets the y component of the current cell

            const scalar y = cellCenter[cellI].component(1);



            // U source term

             const scalar S_u = pi*(-1.6*(2.0*Foam::sin(0.8*pi*x) - 0.6)*Foam::cos(0.8*pi*y) + 1.28*pi*Foam::sin(0.8*pi*y))



             // V source term

           const scalar S_v = -pi*(1.6*(2.0*Foam::sin(0.8*pi*y) + 1.0)*Foam::cos(0.8*pi*x) + 1.28*pi*Foam::sin(0.8*pi*x))



            USource[cellI][0]=-V[cellI]*(S_u);

            USource[cellI][1]=-V[cellI]*(S_v);

            USource[cellI][2]= 0.0;

        };

        

    #};

}

Can anyone confirm this implementation?

Does anyone have experience with this in OpenFOAM? Any advices?

Quote:

Originally Posted by Shibi (Post 802533)

Code:

MomentumSource

{

    type            vectorCodedSource;

    selectionMode   all;

    fields          (U);



    // Name of the coded source

    name            sourceTime;



    codeInclude

    #{



    #};



    codeCorrect

    #{



    #};

       



    codeConstrain

    #{



    #};

    codeAddSup

    #{



        // Info: Include sources into equation



        // Gets the cell volumes of the mesh

        const scalarField& V = mesh_.V();



        // Gets the vector containing cell centers position of the mesh

        const volVectorField& cellCenter = mesh().C();



        // Gets the equation source term

        vectorField& USource = eqn.source();



        const scalar pi = constant::mathematical::pi;



        // Loop over the source term

        forAll(USource, cellI)

        {    

            // Gets the x component of the current cell

            const scalar x = cellCenter[cellI].component(0);



            // Gets the y component of the current cell

            const scalar y = cellCenter[cellI].component(1);



            // U source term

             const scalar S_u = pi*(-1.6*(2.0*Foam::sin(0.8*pi*x) - 0.6)*Foam::cos(0.8*pi*y) + 1.28*pi*Foam::sin(0.8*pi*y))



             // V source term

           const scalar S_v = -pi*(1.6*(2.0*Foam::sin(0.8*pi*y) + 1.0)*Foam::cos(0.8*pi*x) + 1.28*pi*Foam::sin(0.8*pi*x))



            USource[cellI][0]=-V[cellI]*(S_u);

            USource[cellI][1]=-V[cellI]*(S_v);

            USource[cellI][2]= 0.0;

        };

        

    #};

}

Can anyone confirm this implementation?

Does anyone have experience with this in OpenFOAM? Any advices?

Hi, did that work for you?