CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

How to modify icoFoam.C to include source term

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 29, 2012, 00:53
Default How to modify icoFoam.C to include source term
  #1
New Member
 
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 5
erncyc is an unknown quantity at this point
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.
erncyc is offline   Reply With Quote

Old   March 29, 2012, 04:06
Default
  #2
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,067
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
Quote:
Originally Posted by erncyc View Post
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.
try this:
Code:
fvVectorMatrix UEqn
        (
            (1+eta)*fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          
        );

        solve(UEqn == -fvc::grad(p));
nimasam is offline   Reply With Quote

Old   March 29, 2012, 04:32
Default How to modify icoFoam.C to include source term
  #3
New Member
 
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 5
erncyc is an unknown quantity at this point
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.94675e-05, Final residual = 2.71278e-105, No Iterations 22
DILUmyPBiCG: Solving for Uy, Initial residual = 7.4114e-05, Final residual = 5.49146e-105, No Iterations 22
DICPCG: Solving for p, Initial residual = 4.46202e-05, Final residual = 8.4033e-17, No Iterations 97
time step continuity errors : sum local = 7.31098e-19, global = -2.32262e-21, cumulative = 1.11216e-17
DICPCG: Solving for p, Initial residual = 1.64486e-06, Final residual = 9.034e-17, No Iterations 92
time step continuity errors : sum local = 7.20432e-19, global = -1.28127e-21, cumulative = 1.11203e-17
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.94634e-05, Final residual = 2.7106e-105, No Iterations 22
DILUmyPBiCG: Solving for Uy, Initial residual = 7.40934e-05, Final residual = 5.84491e-105, No Iterations 22
DICPCG: Solving for p, Initial residual = 4.46209e-05, Final residual = 8.40861e-17, No Iterations 97
time step continuity errors : sum local = 7.16085e-19, global = 4.24348e-21, cumulative = 1.11246e-17
DICPCG: Solving for p, Initial residual = 1.64556e-06, Final residual = 9.18026e-17, No Iterations 92
time step continuity errors : sum local = 7.14218e-19, global = 3.90334e-21, cumulative = 1.11285e-17
ExecutionTime = 19.96 s ClockTime = 35 s
erncyc is offline   Reply With Quote

Old   March 29, 2012, 05:22
Default
  #4
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,067
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
Quote:
Originally Posted by erncyc View Post
but the result i get is really weird...the solver seems to get 'stuck'
whats wrong with the result? do u mean high iteration numbers?
nimasam is offline   Reply With Quote

Old   March 29, 2012, 05:28
Default
  #5
New Member
 
ernest
Join Date: Jun 2010
Posts: 21
Rep Power: 5
erncyc is an unknown quantity at this point
yes i get very low residuals as soon as the iterations begin and they do not reduce significantly
erncyc is offline   Reply With Quote

Old   March 29, 2012, 05:36
Default
  #6
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,067
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
i have no idea! i dont think it is a cause of error
nimasam is offline   Reply With Quote

Old   February 18, 2013, 09:59
Default
  #7
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 72
Blog Entries: 1
Rep Power: 5
tfuwa is on a distinguished road
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));
where, the Source is defined as,

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) )
    );
This is ok. No error during the compile.

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) )
    );
The compiler complains that:

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_64-linux-gnu/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’ [-Wunused-variable]
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:11:10: warning: unused variable ‘transonic’ [-Wunused-variable]
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:14:9: warning: unused variable ‘nOuterCorr’ [-Wunused-variable]
make: *** [Make/linux64GccDPOpt/icoScalarSFoam.o] Error 1
Could anyone let me know the reason? Thanks.
__________________
Kind regards,

Albert
tfuwa is offline   Reply With Quote

Old   February 18, 2013, 10:10
Default
  #8
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,067
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
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:
Foam::cos (........)
__________________
Training Course on OpenFOAM at (http://www.isme.ir/)
My Weblog (http://openfoam.blogfa.com/)
nimasam is offline   Reply With Quote

Old   February 18, 2013, 10:38
Default
  #9
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 72
Blog Entries: 1
Rep Power: 5
tfuwa is on a distinguished road
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
tfuwa is offline   Reply With Quote

Old   February 18, 2013, 12:38
Default
  #10
Senior Member
 
Nima Sam
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,067
Blog Entries: 1
Rep Power: 14
nimasam is on a distinguished road
Send a message via Yahoo to nimasam
there are two options:
1- define it from Dictionary, for example look at interFoam (dambreak test case) to see how to define a non-uniform 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
__________________
Training Course on OpenFOAM at (http://www.isme.ir/)
My Weblog (http://openfoam.blogfa.com/)
nimasam is offline   Reply With Quote

Old   April 22, 2014, 09:00
Default
  #11
New Member
 
houwy
Join Date: Nov 2013
Posts: 20
Rep Power: 2
houwy is on a distinguished road
Quote:
Originally Posted by tfuwa View Post
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?
Have you solved your problem? I want to do something like yours. I need some help!
houwy is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
momentum source term zwdi FLUENT 13 December 5, 2013 17:58
GPU Linear Solvers for OpenFOAM gocarts OpenFOAM Announcements from Other Sources 35 March 1, 2012 20:41
ATTENTION! Reliability problems in CFX 5.7 Joseph CFX 14 April 20, 2010 15:45
Compiling gmshFoam with OpenFOAM-1.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 11:55


All times are GMT -4. The time now is 05:51.