CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Add soure terms to the momentum equation

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 26, 2021, 16:02
Default Add soure terms to the momentum equation
  #1
Member
 
Join Date: Feb 2020
Posts: 90
Rep Power: 6
Shibi is on a distinguished road
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?
Shibi is offline   Reply With Quote

Old   June 23, 2022, 16:27
Default
  #2
Member
 
Mahmoud
Join Date: Nov 2020
Location: United Kingdom
Posts: 43
Rep Power: 5
Mahmoud Abbaszadeh is on a distinguished road
Quote:
Originally Posted by Shibi View Post
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?
Mahmoud Abbaszadeh is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
How to add a body force to a momentum equation anon_q OpenFOAM Programming & Development 25 May 11, 2021 16:39
Splitting Momentum Eqns and Adding Extra Terms to Each ErenC OpenFOAM 13 May 23, 2020 00:30
Adding two terms in the momentum equation in in chtMultiRegionSimpleFoam zahraa OpenFOAM Pre-Processing 0 June 13, 2015 13:42
Add the constant pressure gradient to the momentum equation tzqfly OpenFOAM Running, Solving & CFD 0 October 6, 2014 05:56
Add source to momentum equation Fan OpenFOAM Programming & Development 0 December 30, 2013 09:35


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