[erncyc]

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.

[nimasam]

Originally Posted by erncyc

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:

fvVectorMatrix UEqn
        (
            (1+eta)*fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
          
        );

        solve(UEqn == -fvc::grad(p));

[erncyc]

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

[nimasam]

Originally Posted by erncyc

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?

[erncyc]

yes i get very low residuals as soon as the iterations begin and they do not reduce significantly

[nimasam]

i have no idea! i dont think it is a cause of error

[tfuwa]

Hi,

I am modifying the icoFoam solver to include a momentum source in y direction.

        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,

        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,

        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:

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.

[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:

Foam::cos (........)

[tfuwa]

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?

[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

[houwy]

Originally Posted by tfuwa

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!

[Luttappy]

Originally Posted by houwy

Have you solved your problem? I want to do something like yours. I need some help!

funkySetFields or setField utility in can be used to initialize this kind of problems.
Please go through Breaking of dam tutorial.