CFD Online Discussion Forums

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)

Shibi April 26, 2021 16:02

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?

Mahmoud Abbaszadeh June 23, 2022 16:27

Quote:

Originally Posted by Shibi (Post 802533)
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?


Hi, did that work for you?


All times are GMT -4. The time now is 18:50.