
[Sponsors] 
March 29, 2012, 00:53 
How to modify icoFoam.C to include source term

#1 
New Member
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 6 
Hi,
I am a novice in the use of OpenFOAM. I was wondering if anyone might be willing to help me modify icoFoam.C so that I can include a source term f, which has the form: f = (eta)*(u_n+1  u_n) / dt u_n+1 and u_n are the current and previous velocities. eta is a variable defined in setFields and created in create.H. I tried the code below but I don't seem to get good results. Indeed if i choose smaller time steps, the solver crashes. fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U)  fvm::laplacian(nu, U) ); solve(UEqn == fvc::grad(p)eta*(fvm::ddt(U))); Please help, thank you. 

March 29, 2012, 04:06 

#2  
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,168
Blog Entries: 1
Rep Power: 15 
Quote:
Code:
fvVectorMatrix UEqn ( (1+eta)*fvm::ddt(U) + fvm::div(phi, U)  fvm::laplacian(nu, U) ); solve(UEqn == fvc::grad(p)); 

March 29, 2012, 04:32 
How to modify icoFoam.C to include source term

#3 
New Member
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 6 
Hi Nima,
Thank you for your prompt reply. I have tried that option but the result i get is really weird...the solver seems to get 'stuck' as can be seen below Courant Number mean: 0.00944072 max: 0.0228351 DILUmyPBiCG: Solving for Ux, Initial residual = 4.94675e05, Final residual = 2.71278e105, No Iterations 22 DILUmyPBiCG: Solving for Uy, Initial residual = 7.4114e05, Final residual = 5.49146e105, No Iterations 22 DICPCG: Solving for p, Initial residual = 4.46202e05, Final residual = 8.4033e17, No Iterations 97 time step continuity errors : sum local = 7.31098e19, global = 2.32262e21, cumulative = 1.11216e17 DICPCG: Solving for p, Initial residual = 1.64486e06, Final residual = 9.034e17, No Iterations 92 time step continuity errors : sum local = 7.20432e19, global = 1.28127e21, cumulative = 1.11203e17 ExecutionTime = 19.94 s ClockTime = 35 s Time = 0.6035 Courant Number mean: 0.00944072 max: 0.0228347 DILUmyPBiCG: Solving for Ux, Initial residual = 4.94634e05, Final residual = 2.7106e105, No Iterations 22 DILUmyPBiCG: Solving for Uy, Initial residual = 7.40934e05, Final residual = 5.84491e105, No Iterations 22 DICPCG: Solving for p, Initial residual = 4.46209e05, Final residual = 8.40861e17, No Iterations 97 time step continuity errors : sum local = 7.16085e19, global = 4.24348e21, cumulative = 1.11246e17 DICPCG: Solving for p, Initial residual = 1.64556e06, Final residual = 9.18026e17, No Iterations 92 time step continuity errors : sum local = 7.14218e19, global = 3.90334e21, cumulative = 1.11285e17 ExecutionTime = 19.96 s ClockTime = 35 s 

March 29, 2012, 05:28 

#5 
New Member
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 6 
yes i get very low residuals as soon as the iterations begin and they do not reduce significantly


February 18, 2013, 10:59 

#7 
Member
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7 
Hi,
I am modifying the icoFoam solver to include a momentum source in y direction. Code:
fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U)  fvm::laplacian(nu, U) + Source ); solve(UEqn == fvc::grad(p)); Code:
volVectorField Source ( IOobject ( "Source", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedVector("Source", dimensionSet(0,1,2,0,0,0,0),vector(0, 100, 0) ) ); However, when I tried to make the Source time dependent as follows, Code:
volVectorField Source ( IOobject ( "Source", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedVector("Source", dimensionSet(0,1,2,0,0,0,0),vector(0, cos(runTime.time().value()), 0) ) ); Code:
createFields.H: In function ‘int main(int, char**)’: createFields.H:62:102: error: call of overloaded ‘cos(const double&)’ is ambiguous createFields.H:62:102: note: candidates are: /usr/include/x86_64linuxgnu/bits/mathcalls.h:64:1: note: double cos(double) /opt/openfoam171/src/OpenFOAM/lnInclude/dimensionedScalar.H:76:19: note: Foam::dimensionedScalar Foam::cos(const dimensionedScalar&) /opt/openfoam171/src/OpenFOAM/lnInclude/Scalar.H:238:1: note: Foam::doubleScalar Foam::cos(Foam::doubleScalar) /opt/openfoam171/src/OpenFOAM/lnInclude/Scalar.H:238:1: note: Foam::floatScalar Foam::cos(Foam::floatScalar) /opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:8:10: warning: unused variable ‘momentumPredictor’ [Wunusedvariable] /opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:11:10: warning: unused variable ‘transonic’ [Wunusedvariable] /opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:14:9: warning: unused variable ‘nOuterCorr’ [Wunusedvariable] make: *** [Make/linux64GccDPOpt/icoScalarSFoam.o] Error 1
__________________
Kind regards, Albert 

February 18, 2013, 11:10 

#8  
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,168
Blog Entries: 1
Rep Power: 15 
as function "cos" is defined in math.h too!, so you should define which one, you want to use, so it would be solve with following structure:
Quote:
__________________
Telegram channel (https://telegram.me/openfoam4Iranian) My Weblog (http://openfoam.blogfa.com/) Training Course on OpenFOAM at (http://www.isme.ir/) 

February 18, 2013, 11:38 

#9 
Member
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 7 
Hi Nima,
Many thanks for your quick reply. The problem is solved. Now, I would like to see the Source term only exists at a specific area. if (x y z) in zone (x1 to x2, y1 to y2, z1 to z2) Source = (0, a, 0) else Source = (0, 0, 0) endif Can you please enlighten me on how to realize this?
__________________
Kind regards, Albert 

February 18, 2013, 13:38 

#10 
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,168
Blog Entries: 1
Rep Power: 15 
there are two options:
1 define it from Dictionary, for example look at interFoam (dambreak test case) to see how to define a nonuniform internalField 2 you can write a pieces of code, you can use mesh.C() which returns position vector then you can check it whether your cell is set up in specific area or not
__________________
Telegram channel (https://telegram.me/openfoam4Iranian) My Weblog (http://openfoam.blogfa.com/) Training Course on OpenFOAM at (http://www.isme.ir/) 

April 22, 2014, 09:00 

#11  
New Member
houwy
Join Date: Nov 2013
Posts: 21
Rep Power: 4 
Quote:


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
momentum source term  zwdi  FLUENT  13  December 5, 2013 18:58 
GPU Linear Solvers for OpenFOAM  gocarts  OpenFOAM Announcements from Other Sources  35  March 1, 2012 21:41 
ATTENTION! Reliability problems in CFX 5.7  Joseph  CFX  14  April 20, 2010 15:45 
Compiling gmshFoam with OpenFOAM1.5  BlGene  Open Source Meshers: Gmsh, Netgen, CGNS, ...  10  August 6, 2009 04:26 
fluent add additional zones for the mesh file  SSL  FLUENT  2  January 26, 2008 12:55 