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

vectorCodedSource with V momentum equation.

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 27, 2019, 08:22
Default vectorCodedSource with V momentum equation.
  #1
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
Hi,


I was trying to modify the V momentum equation with a source term. I tried it with the vectorCodedSource term specification as given in the below code snippet. The issue is that the case runs but the source terms are not being included while solving the problem. i.e the solution I get is similar to the solution I get without employing the source term.

I am using simpleFoam (openFOAM v1906) as my solver and the test case I am solving corresponds to a lid driven cavity with manufactured solutions. The reference paper is :
Shih, T. M., Tan, C. H., & Hwang, B. C. (1989). Effects of grid staggering on numerical schemes. International Journal for Numerical Methods in Fluids, 9(2), 193212. https://doi.org/10.1002/fld.1650090206



Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system" ;
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


MomentumSource 

{ type vectorCodedSource ;
  name MomentumSource ;
  active on ;  
  fields   (U) ;
  selectionMode all ;
  redirectType velocitySource; 
   
  codeInclude
  #{

  #};

  codeCorrect
  #{

    Info<< "CodeCorrect" << endl << "\n" ; 
  
   #};

  codeConstrain
  #{
   
    Info<< "CodeConstrain" << endl << "\n" ; 

  #};


  codeAddSup
  #{ 
     
     Info<< "codeAddSup" << endl << "\n" ; 

     const scalarField&    V  = mesh_.V() ;  // Cell volume  

     vectorField& VSource = eqn.source() ;   // Equation source definition

     const vectorField& cell  = mesh_.C()  ;
     const scalarField& cellx = cell.component(0) ; 
     const scalarField& celly = cell.component(1) ; 
     

     int n = cell.size() ;  
     double arg[n] = {} ; 
     double v1[n] = {} ;   
     double v2[n] = {} ;  
     double v3[n] = {} ;  
     double v4[n] = {} ;  
     double v5[n] = {} ;  
     double v6[n] = {} ;  
     double v7[n] = {} ;  
     double v8[n] = {} ;  
     double v9[n] = {} ;  
     double v10[n] = {} ;  

Info<<"VSource"<<VSource<<"\n" ; 

     forAll(cell, i) 
     {

     v1[i] = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4) - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;


     v2[i] = -768*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 8*pow(cellx[i],6) - 10*pow(cellx[i],5) + 8*pow(cellx[i],4) - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],5) ; 


     v3[i] = 38.4*pow(cellx[i],5) ;  

     v4[i] = -96*pow(cellx[i],4) ; 

     v5[i] = 256*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 8*pow(cellx[i],6) - 10*pow(cellx[i],5) + 8*pow(cellx[i],4) - 4*pow(cellx[i],3) +  pow(cellx[i],2) ) * pow(celly[i],3) ;

     v6[i] = -64*pow(cellx[i],3) ; 


     v7[i] =  96*(8*pow(cellx[i],3) - 12*pow(cellx[i],2) + 2*pow(cellx[i],1) + 1)*pow(celly[i],2) ;

     v8[i] = 192*pow(cellx[i],2) ;

     v9[i] = -128*(pow(cellx[i],8) - 4*pow(cellx[i],7) + 6*pow(cellx[i],6) - 4*pow(cellx[i],5) + pow(cellx[i],4))*pow(celly[i],1);

     v10[i] = -64*cellx[i] ; 
 
      arg[i] = (v1[i] + v2[i] + v3[i] + v4[i] + v5[i] + v6[i] + v7[i] + v8[i] + v9[i] + v10[i])*V[i]  ; 

      VSource[i] = vector(0,arg[i],0) ;

     } 



  #} ; 


code
        #{
            $codeInclude
            $codeCorrect
            $codeConstrain
            $codeAddSup
           $codeSetValue
        #};
  }
Any ideas on tackling this issue ?

Thanks in advance,
Praphul
praphul is offline   Reply With Quote

Old   August 28, 2019, 06:57
Default
  #2
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
I tried the same with the following command. Just to see whether it invokes the source term but it doesnt do so. Am I doing anything wrong?

Code:
 vectorField& Usource = eqn.source();

       Usource += vector(0, 10000, 0);
praphul is offline   Reply With Quote

Old   August 29, 2019, 21:14
Default
  #3
Senior Member
 
Pete Bachant
Join Date: Jun 2012
Location: Boston, MA
Posts: 171
Rep Power: 11
pbachant is on a distinguished road
You may want to try creating a field an adding it to eqn directly, i.e., eqn += myField, or you could try adding to the equation's source rather than setting the values explicitly. I don't know enough about the inner workings to provide more insight than that!
__________________
Home | Twitter | GitHub
pbachant is offline   Reply With Quote

Old   September 1, 2019, 03:21
Default
  #4
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
Thank you pbachant. I tried setting the field directly to the equation as you had mentioned. Unfortunately, It does not work. But, I Thank you for your time and consideration.
praphul is offline   Reply With Quote

Old   September 1, 2019, 12:59
Default
  #5
Member
 
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6
ELwardi is on a distinguished road
Well, I see you tagged the thread with openfoamv2019; I don't have OF installed (Just fragments of foam-extend 4 actually ). So I'll try to help with no promise that the code will even compile:


1. Avoid using C ARRAYS. At least, use Foam::List if absolutely required .


2. You create way too much temp vars; Try getting to the point faster: You want to modify the second component in the source of UEqn? get to it immediately in the forAll loop.


3. Always ADD values to source, don't assign to it !!!


It has been a while since I last medeled with fvOptions, but I will share what I would write to accomplish this (Of course, I welcome any improvements/suggestions):
Code:
codedSource
{
    type            coded;
    selectionMode   all;

    fields          (U);
    name            MomentumSource;

    codeAddSup
    #{
        const scalarField& V = mesh_.V();
        scalarField& USource = eqn.source();
        const scalarField& cellx = mesh_.C().component(0) ; 
        const scalarField& celly = mesh_.C().component(0) ;
        
        forAll(USource, i)
        {
            scalar v1 = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4)
                             - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;
            // .... and the other 9 scalars
            
            // Please ADD to the source, don't overwrite it.
            // Overwriting it means you don't care about source terms in 
            // DESCRETISATION OPERATORS (laplacian, div ...)
            USource[i][1] += (v1+v2+...)*V[i];
        }
    #};
}
I hope it fixes your problem. Maybe you didn't notice the change by overwriting the source.
ELwardi is offline   Reply With Quote

Old   September 3, 2019, 08:25
Default
  #6
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
Elwardi, I tried ur suggestion and unfortunately that doesn't work too. It seems I am stuck with this problem for eternity.

1) I am curious to know on ur suggestion as to why I should not use C arrays and instead I should use list.



2) The Velocity field is supposed to be a vector quantity. Therefore, I am supposed to use vectorCodedSource isnt it ? I guess the way you have specified it might be because you use a different version of OpenFoam.


3) Why is it said that we should always add to the source rather than assigning it ?

I am attaching the zip file of my problem setup with this mail. May be if you or someone else can look into that, it can be helpful .

Regards,
Praphul
Attached Files
File Type: zip LidDrivenCavity.zip (29.5 KB, 6 views)
praphul is offline   Reply With Quote

Old   September 3, 2019, 11:25
Default
  #7
Member
 
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6
ELwardi is on a distinguished road
Ok, let me clarify things a bit:


1. It's not like C arrays are particularly bad for this situation, but it's good to jump out of the habit of using them (No bounds checking, possibility to overflow the stack ...) http://www.cplusplus.com/forum/articles/17108/


2. I meant "vectorField& USource = eqn.source();", that was a copy-paste typo.


3. I was considering that equation source may be already be populated by the time fvOptions kicks in!
Suppose your fvm::laplacian generates an explicit term (fvm::ddt certainly does in transient simulations). That term will be stored in the source vector. By the time your additional source is applied, if you overwrite what was there with your source's value, you're no longer solving the problem you want to solve
But it appears recent fvOptions generate their own hole "eqn" then adds it the main equation. I almost always access sources directly so i got confused


4. I had problems compiling your U boundary conditions (writeEntry seems to be dropped in OF7) but I managed to get it working (At least there is a change) with basic U BCs. Let me know if the problem persists for you
Attached Files
File Type: zip LidDrivenCavity.zip (31.6 KB, 7 views)
ELwardi is offline   Reply With Quote

Old   September 4, 2019, 03:28
Default
  #8
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
Thank you for clarifying my doubts Elwardi.



1) However, in the file 0/U you had attached I see that you had specified zeroGradient boundary condition for the walls which I think should be changed to noSlip condition.

2) I do not understand why the celly component was specified like this is in the fvOptions.
Code:
        const scalarField& celly = mesh_.C().component(0) ;
y component should have an index 1 isnt it ? i.e
Code:
        const scalarField& celly = mesh_.C().component(1) ;
I changed both these criterion and ran the case. I found that the solution I get is exactly the same as with fvOptions and without fvOptions. I dont understand why that is so. I have attached graphs for the same.
Attached Images
File Type: jpeg With FvOptions.jpeg (33.7 KB, 33 views)
File Type: jpeg Without FvOptions.jpeg (34.4 KB, 13 views)
praphul is offline   Reply With Quote

Old   September 4, 2019, 08:47
Default
  #9
Member
 
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6
ELwardi is on a distinguished road
After fixing the issues you mentioned, I got slower convergence with the source enabled (340 iterations vs 330 iters in a x4 finer mesh). Maybe the source is too week to notice the change, try to multiply USource with 10 in all cells, there will be a clear change to notice. Revise your vi quantities, there might be a bug in there!
ELwardi is offline   Reply With Quote

Old   September 7, 2019, 03:01
Default
  #10
New Member
 
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9
praphul is on a distinguished road
In the particular case of manufactured solution which I solve, I assume a solution for u, v and plug it in to the governing equations to obtain the source terms for the respective equation. In this case, I get the x momentum source term to be zero and the y momentum source term to be what I have specified in the fvOptions. Then I solve the equations again to get back the solution I have assumed. Its like a merry-go-round but it is used to verify the code.

To get the same solution as the assumed one, the boundary conditions needs to be in a format I had specified in the zip file earlier.

I dont think the source term comes out to be small as I had tried printing out the value of sources using Info command and it happen to be finite value. The research paper suggests that the source term at x = 0.5 and y = 0.5 and Reynolds number of 1 equals 3.356250 and I get the same value. In this example, this value equals source term divded by the volume of the cell. I am supposed to get back the original solution I had assumed unless there is some problem with the code.

I had cross checked the vi terms too. Infact I had missed a term i.e 96*(2x-1)*y^4 in the fvOption file I had specified, I corrected it but the solution remains same. Still wondering what the issue could be.
praphul is offline   Reply With Quote

Old   September 7, 2019, 06:52
Default
  #11
Member
 
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6
ELwardi is on a distinguished road
1. Double check your units! I think UEqn has [0 4 -2 0 0 0 0]



2. You can use:
debugSwitches{
vectorCodedSource 1;

}
in controlDict to verify when ::addSup method is called and whether it does what you want when you want.
ELwardi is offline   Reply With Quote

Old   September 8, 2019, 07:45
Default
  #12
Member
 
CFD USER
Join Date: May 2019
Posts: 40
Rep Power: 4
CFD_10 is on a distinguished road
Quote:
Originally Posted by ELwardi View Post
1. Double check your units! I think UEqn has [0 4 -2 0 0 0 0]



2. You can use:
debugSwitches{
vectorCodedSource 1;

}
in controlDict to verify when ::addSup method is called and whether it does what you want when you want.
I cannot see vectorCodedSource in the global controlDict of (OF5, OF6, and OF 1906).
Code:
grep -i 'codedVectorSource' $FOAM_SRC/controlDict
This will return nothing.
CFD_10 is offline   Reply With Quote

Old   September 8, 2019, 07:50
Default
  #13
Member
 
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6
ELwardi is on a distinguished road
`simpleFoam -listSwitches | grep Source`
Tells me there is one in OF7; there should be something similar eg. codedSource, run the command and see if you can find any useful switch
ELwardi is offline   Reply With Quote

Old   September 8, 2019, 08:24
Default
  #14
Member
 
CFD USER
Join Date: May 2019
Posts: 40
Rep Power: 4
CFD_10 is on a distinguished road
Quote:
Originally Posted by ELwardi View Post
`simpleFoam -listSwitches | grep Source`
Tells me there is one in OF7; there should be something similar eg. codedSource, run the command and see if you can find any useful switch
According to the solver help, the option -listSwitches : List switches declared in libraries but not set in etc/controlDict
CFD_10 is offline   Reply With Quote

Reply

Tags
openfoamv1906, programming, vectorcodedsource

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
add a pressure drop term in the momentum equation a.lone FLUENT 0 July 3, 2019 06:48
Domain Reference Pressure and mass flow inlet boundary AdidaKK CFX 75 August 20, 2018 05:37
mass flow in is not equal to mass flow out saii CFX 12 March 19, 2018 05:21
Coefficients discretized momentum equation michujo Main CFD Forum 4 June 20, 2012 01:33
Discretization of momentum equation query siw CFX 0 June 20, 2011 08:38


All times are GMT -4. The time now is 14:06.