
[Sponsors] 
How to modify icoFoam.C to include source term 

LinkBack  Thread Tools  Search this Thread  Display Modes 
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: 15 
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,266
Blog Entries: 1
Rep Power: 24 
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: 15 
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: 15 
yes i get very low residuals as soon as the iterations begin and they do not reduce significantly


February 18, 2013, 09:59 

#7 
Member
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 15 
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, 10:10 

#8  
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,266
Blog Entries: 1
Rep Power: 24 
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:
__________________
My Personal Website (http://nimasamkhaniani.ir/) Telegram channel (https://t.me/cfd_foam) 

February 18, 2013, 10:38 

#9 
Member
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 76
Blog Entries: 1
Rep Power: 15 
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, 12:38 

#10 
Senior Member
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,266
Blog Entries: 1
Rep Power: 24 
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
__________________
My Personal Website (http://nimasamkhaniani.ir/) Telegram channel (https://t.me/cfd_foam) 

April 22, 2014, 09:00 

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


June 10, 2018, 01:38 
setField or funkySetField

#12  
Member
HK
Join Date: Oct 2015
Location: Madras
Posts: 31
Rep Power: 10 
Quote:
Please go through Breaking of dam tutorial. 

Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
GPU Linear Solvers for OpenFOAM  gocarts  OpenFOAM Announcements from Other Sources  37  August 17, 2022 14:22 
momentum source term  zwdi  FLUENT  14  June 27, 2017 15:40 
ATTENTION! Reliability problems in CFX 5.7  Joseph  CFX  14  April 20, 2010 15:45 
[Gmsh] Compiling gmshFoam with OpenFOAM1.5  BlGene  OpenFOAM Meshing & Mesh Conversion  10  August 6, 2009 04:26 
fluent add additional zones for the mesh file  SSL  FLUENT  2  January 26, 2008 11:55 