
[Sponsors] 
August 27, 2019, 08:22 
vectorCodedSource with V momentum equation.

#1 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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), 193–212. 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 #}; } Thanks in advance, Praphul 

August 28, 2019, 06:57 

#2 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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); 

August 29, 2019, 21:14 

#3 
Senior Member
Pete Bachant
Join Date: Jun 2012
Location: Boston, MA
Posts: 171
Rep Power: 11 
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!


September 1, 2019, 03:21 

#4 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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.


September 1, 2019, 12:59 

#5 
Member
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6 
Well, I see you tagged the thread with openfoamv2019; I don't have OF installed (Just fragments of foamextend 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]; } #}; } 

September 3, 2019, 08:25 

#6 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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 

September 3, 2019, 11:25 

#7 
Member
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6 
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 copypaste 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 

September 4, 2019, 03:28 

#8 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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) ; Code:
const scalarField& celly = mesh_.C().component(1) ; 

September 4, 2019, 08:47 

#9 
Member
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6 
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!


September 7, 2019, 03:01 

#10 
New Member
Praphul T
Join Date: Dec 2013
Location: Bangalore, India
Posts: 23
Rep Power: 9 
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 merrygoround 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*(2x1)*y^4 in the fvOption file I had specified, I corrected it but the solution remains same. Still wondering what the issue could be. 

September 7, 2019, 06:52 

#11 
Member
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6 
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. 

September 8, 2019, 07:45 

#12  
Member
CFD USER
Join Date: May 2019
Posts: 40
Rep Power: 4 
Quote:
Code:
grep i 'codedVectorSource' $FOAM_SRC/controlDict 

September 8, 2019, 07:50 

#13 
Member
Elwardi Fadeli
Join Date: Dec 2016
Location: Boumerdes, Algeria
Posts: 39
Rep Power: 6 
`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 

September 8, 2019, 08:24 

#14 
Member
CFD USER
Join Date: May 2019
Posts: 40
Rep Power: 4 
According to the solver help, the option listSwitches : List switches declared in libraries but not set in etc/controlDict


Tags 
openfoamv1906, programming, vectorcodedsource 
Thread Tools  Search this Thread 
Display Modes  


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 