CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Adding source term with fvOptions to k-e model

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 21, 2020, 16:34
Default Adding source term with fvOptions to k-e model
  #1
Member
 
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2
Tibo99 is on a distinguished road
Hi everyone!

I got some issue to implement a source term via fvOptions to the k-e model.

The source term is apply on epsilon(only) to the k-e turbulence model for 1D (x=1mm, y=150mm, z=1mm) ABL simulation in order to simplify the simulation.

Y+ is around 50. The y component I used in the formula is the center of the cells so, y-coordinate of internal cells (distance from the wall to the center of each cells,see text file, 3rd column).

From the results, the calculation by the equation "epsilonSource[i]" are good since I printed out the results in a text file and verified the validity but the source term doesn't seems to be applied correctly to the domain(see picture). The boundary wall seems to work (in the epsilon patch I choose the BC "epsilonWallFunction"). The simulation stop after 11 step of iteration...

If someone has a clue, an idea what could be the issue, it would be appreciated.

Do you guys think the source term is applied automatically in the epsilon-BC as well? If don't, should I?

Best Regards,

Paper where I found the term source: https://www.sciencedirect.com/scienc...6761051100002X
https://www.researchgate.net/publica...n_of_ABL_flows

fvOptions file:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1906                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

epsilonSource
{
    type            scalarCodedSource;
    selectionMode   all;
    active			true;

    fields          (epsilon);
    name			epsilonSource;
    
    
    codeInclude
    #{
		#include "OFstream.H"
	#};
    
    codeAddSup
    #{
	const Time& time = mesh_.time();
				
        const scalar us_ = 0.76;
        const scalar D1k_ = 1.44;
        const scalar D2k_ = 1.92;
        const scalar kappa_ = 0.41;
        const scalar sigmaEps_ = 1.3;
        const scalar Cmu_ = 0.09;
        const scalar z0_ = 3.3474e-05;
        
        const vectorField& CellC = mesh_.C(); 
        
		
	scalarField& epsilonSource = eqn.source();
				
        // Apply the source
        forAll(CellC, i)
        {
          							
             
               epsilonSource[i] += (pow(us_,4)/(pow((CellC[i].y()+z0_),2)))
               *(((D2k_-D1k_)*sqrt(Cmu_)/pow(kappa_,2))-(1/sigmaEps_));
            
          		   			
		std::ofstream file;
		file.open ("Epsilonterm.txt", std::ofstream::out | std::ofstream::app);
		file << time.value() << "	" << epsilonSource[i] << "	" << CellC[i].y() << std::endl << "\n";
		file.close();				
        
        };
        
        
    #};
    
    codeCorrect
    #{
    #};
    
    codeConstrain
    #{
    #};
}
// ************************************************************************* //
Log file:
Code:
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1906                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : v1906 OPENFOAM=1906
Arch   : "LSB;label=32;scalar=64"
Exec   : simpleFoam
Date   : May 21 2020
Time   : 17:23:59
Host   : rene-OptiPlex-790
PID    : 16903
I/O    : uncollated
Case   : /home/rene/OpenFOAM/rene-v1906/run/1D_E1
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: convergence criteria
    field U	 tolerance 1e-06
    field k	 tolerance 1e-06
    field epsilon	 tolerance 1e-06
    field omega	 tolerance 1e-06
    field nuTilda	 tolerance 1e-06
    field T	 tolerance 1e-06
    field p_rgh	 tolerance 1e-06
    field p	 tolerance 1e-06

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    label           "Standard high-Re k-\u03B5";
    fieldMaps
    {
        k               k;
        epsilon         epsilon;
        nut             nut;
    }
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    E               9.8;
    sigmak          1;
    sigmaEps        1.3;
    C3              0;
}

No MRF models present

Creating finite volume options from "system/fvOptions"

Selecting finite volume options type scalarCodedSource
    Source: epsilonSource
    - selecting all cells
    - selected 150 cell(s) with volume 1.5e-07

Starting time loop

turbulenceFields turbulenceFields1: storing fields:
    turbulenceProperties:R
    turbulenceProperties:I
    turbulenceProperties:L

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 5.781665123e-06, Final residual = 2.320453622e-10, No Iterations 8
smoothSolver:  Solving for Uy, Initial residual = 0.9089618075, Final residual = 5.582492529e-05, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 1.786751413e-15, Final residual = 4.823728453e-16, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.7530448709, Final residual = 5.369605972e-05, No Iterations 21
time step continuity errors : sum local = 3.673528501e-19, global = 3.673528501e-19, cumulative = 3.673528501e-19
Using dynamicCode for fvOption::epsilonSource at line -1 in "/home/rene/OpenFOAM/rene-v1906/run/1D_E1/system/fvOptions.epsilonSource"
Creating new library in "dynamicCode/epsilonSource/platforms/linux64GccDPInt32Opt/lib/libepsilonSource_273d0b6b78f639d7c93a407718428cb1cf5e4e35.so"
Invoking wmake libso /home/rene/OpenFOAM/rene-v1906/run/1D_E1/dynamicCode/epsilonSource
wmake libso /home/rene/OpenFOAM/rene-v1906/run/1D_E1/dynamicCode/epsilonSource
    dep: codedFvOptionTemplate.C
    Ctoo: codedFvOptionTemplate.C
    ld: /home/rene/OpenFOAM/rene-v1906/run/1D_E1/dynamicCode/epsilonSource/../platforms/linux64GccDPInt32Opt/lib/libepsilonSource_273d0b6b78f639d7c93a407718428cb1cf5e4e35.so
Selecting finite volume options type epsilonSource
    Source: epsilonSource
    - selecting all cells
    - selected 150 cell(s) with volume 1.5e-07
smoothSolver:  Solving for epsilon, Initial residual = 0.9998814658, Final residual = 2.827019512e-05, No Iterations 8
bounding epsilon, min: -506399185.4 max: 2140.790268 average: -5800585.107
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 9.257761502e-05, No Iterations 7
ExecutionTime = 0.06 s  ClockTime = 3 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.0002397693829, Final residual = 8.113253956e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.01029908582, Final residual = 4.358351871e-07, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.0629515603, Final residual = 2.569515862e-06, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.7143023531, Final residual = 104040663.2, No Iterations 1000
time step continuity errors : sum local = 1.213998486e-05, global = 1.400000834e-09, cumulative = 1.400000834e-09
smoothSolver:  Solving for epsilon, Initial residual = 0.003063504061, Final residual = 2.350539837e-07, No Iterations 6
bounding epsilon, min: -4.792808503e-11 max: 1947.869938 average: 13.85865537
smoothSolver:  Solving for k, Initial residual = 0.397705588, Final residual = 2.981436582e-05, No Iterations 7
ExecutionTime = 0.3 s  ClockTime = 3 s

Time = 3

smoothSolver:  Solving for Ux, Initial residual = 6.542876533e-05, Final residual = 3.414747811e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.008345991815, Final residual = 6.916702663e-07, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.4207610645, Final residual = 1.767895306e-05, No Iterations 7
GAMG:  Solving for p, Initial residual = 0.9875079438, Final residual = 2.876784968e-05, No Iterations 5
time step continuity errors : sum local = 72.2530223, global = 7.186872358e-10, cumulative = 2.11868807e-09
smoothSolver:  Solving for epsilon, Initial residual = 1.779790796e-10, Final residual = 2.332839125e-11, No Iterations 1
bounding epsilon, min: -0.008365059191 max: 926008400.3 average: 10967512.72
smoothSolver:  Solving for k, Initial residual = 0.2500923335, Final residual = 2.109144553e-05, No Iterations 7
ExecutionTime = 0.3 s  ClockTime = 3 s

Time = 4

smoothSolver:  Solving for Ux, Initial residual = 0.0001047864544, Final residual = 3.802370649e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 1.067786309e-09, Final residual = 9.803739005e-11, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 5.08984362e-07, Final residual = 3.922613063e-11, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.9999368938, Final residual = 9.488941883e-05, No Iterations 974
time step continuity errors : sum local = 94.93563983, global = 2.048071334e-05, cumulative = 2.048283203e-05
smoothSolver:  Solving for epsilon, Initial residual = 1.464161316e-10, Final residual = 2.944491158e-11, No Iterations 1
bounding epsilon, min: -2.276798111e-11 max: 1591232665 average: 14897422.52
smoothSolver:  Solving for k, Initial residual = 4.594718889e-09, Final residual = 8.271137915e-11, No Iterations 3
ExecutionTime = 0.36 s  ClockTime = 3 s

Time = 5

smoothSolver:  Solving for Ux, Initial residual = 5.71021682e-05, Final residual = 3.086792287e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 1.395646902e-09, Final residual = 6.012380828e-11, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 1.899201834e-10, Final residual = 3.854897013e-11, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.9998158918, Final residual = 6.610612207e-05, No Iterations 10
time step continuity errors : sum local = 15182.47696, global = -0.005463523176, cumulative = -0.005443040344
smoothSolver:  Solving for epsilon, Initial residual = 5.515123834e-09, Final residual = 3.612469016e-11, No Iterations 4
bounding epsilon, min: -1.555379741e-11 max: 3.206587777e+12 average: 2.633791325e+10
smoothSolver:  Solving for k, Initial residual = 8.865456493e-07, Final residual = 6.354148454e-11, No Iterations 7
ExecutionTime = 0.37 s  ClockTime = 3 s

Time = 6

smoothSolver:  Solving for Ux, Initial residual = 4.903075674e-05, Final residual = 2.76209214e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.001152352045, Final residual = 3.81812229e-08, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.00082024583, Final residual = 5.758955053e-08, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.9912377818, Final residual = 6.390959052e-05, No Iterations 10
time step continuity errors : sum local = 5758.516424, global = 183.0669591, cumulative = 183.0615161
smoothSolver:  Solving for epsilon, Initial residual = 0.0003505345651, Final residual = 1.672863807e-08, No Iterations 7
bounding epsilon, min: -1.131701388e-11 max: 7.930001666e+12 average: 8.597395232e+10
smoothSolver:  Solving for k, Initial residual = 0.01084036878, Final residual = 7.70259312e-07, No Iterations 7
ExecutionTime = 0.38 s  ClockTime = 3 s

Time = 7

smoothSolver:  Solving for Ux, Initial residual = 0.0001017155284, Final residual = 4.623390674e-09, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.005494247008, Final residual = 1.851745399e-07, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.01430007137, Final residual = 4.569324173e-07, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.8459483358, Final residual = 3.237544819e-05, No Iterations 12
time step continuity errors : sum local = 6614.508604, global = 13.01227121, cumulative = 196.0737873
smoothSolver:  Solving for epsilon, Initial residual = 7.049825885e-06, Final residual = 3.746246606e-10, No Iterations 7
bounding epsilon, min: -7.098149768e-12 max: 1.287364284e+14 average: 9.485271478e+11
smoothSolver:  Solving for k, Initial residual = 0.001794322734, Final residual = 1.348391127e-07, No Iterations 7
ExecutionTime = 0.38 s  ClockTime = 3 s

Time = 8

smoothSolver:  Solving for Ux, Initial residual = 0.001637598567, Final residual = 6.434227911e-08, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.1603902233, Final residual = 4.344377787e-06, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 0.0002409532143, Final residual = 1.363044774e-08, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.623732765, Final residual = 14.43795073, No Iterations 1000
time step continuity errors : sum local = 1923661095, global = 32891581.74, cumulative = 32891777.81
smoothSolver:  Solving for epsilon, Initial residual = 8.42542935e-11, Final residual = 2.006210096e-11, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.4804684285, Final residual = 4.048431869e-05, No Iterations 7
ExecutionTime = 0.59 s  ClockTime = 4 s

Time = 9

smoothSolver:  Solving for Ux, Initial residual = 0.000961475335, Final residual = 4.019731076e-08, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.05898306748, Final residual = 1.477984489e-06, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 3.178084763e-08, Final residual = 2.804660479e-11, No Iterations 4
GAMG:  Solving for p, Initial residual = 0.9148160798, Final residual = 0.261447627, No Iterations 1000
time step continuity errors : sum local = 7840704.186, global = 7825413.287, cumulative = 40717191.1
smoothSolver:  Solving for epsilon, Initial residual = 3.413948522e-17, Final residual = 7.529783842e-18, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.2539284515, Final residual = 2.201685003e-05, No Iterations 7
ExecutionTime = 0.84 s  ClockTime = 4 s

Time = 10

smoothSolver:  Solving for Ux, Initial residual = 0.0009077776947, Final residual = 3.482148807e-08, No Iterations 7
smoothSolver:  Solving for Uy, Initial residual = 0.1731602652, Final residual = 7.857890645e-06, No Iterations 7
smoothSolver:  Solving for Uz, Initial residual = 2.73445177e-08, Final residual = 2.887628211e-11, No Iterations 4
GAMG:  Solving for p, Initial residual = 0.000899333075, Final residual = 2.607603027e+68, No Iterations 1000
time step continuity errors : sum local = 9.804016467e+79, global = 4.47494892e+74, cumulative = 4.47494892e+74
smoothSolver:  Solving for epsilon, Initial residual = 1, Final residual = 2.155609824e-11, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 9.253973429e-05, No Iterations 7
ExecutionTime = 1.33 s  ClockTime = 4 s

Time = 11

smoothSolver:  Solving for Ux, Initial residual = 1.40996752e-13, Final residual = 5.783231539e-15, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 9.056049294e-07, Final residual = 4.438814615e-11, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 1.40996752e-13, Final residual = 5.783231539e-15, No Iterations 1
Attached Images
File Type: jpg Epsilon.jpg (97.1 KB, 8 views)
File Type: jpg nut.jpg (101.1 KB, 6 views)
File Type: jpg Ux.jpg (96.5 KB, 5 views)
File Type: png Esource.png (3.3 KB, 5 views)
Attached Files
File Type: txt Epsilonterm.txt (26.3 KB, 4 views)

Last edited by Tibo99; May 24, 2020 at 10:22.
Tibo99 is offline   Reply With Quote

Old   May 28, 2020, 07:09
Default
  #2
Member
 
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2
superkelle is on a distinguished road
Hi, do you get an floating point error after t=11? I am not sure about the "+=" maybe try simple "=".
superkelle is offline   Reply With Quote

Old   May 28, 2020, 12:49
Default
  #3
Member
 
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2
Tibo99 is on a distinguished road
I already tried ''+='', ''='' and ''-='', and I got the same error.

That is the last part of the log:

Code:
Time = 11

smoothSolver:  Solving for Ux, Initial residual = 1.40996752e-13, Final residual = 5.783231539e-15, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 9.056049294e-07, Final residual = 4.438814615e-11, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 1.40996752e-13, Final residual = 5.783231539e-15, No Iterations 1
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::scalarProduct<double, double>::type Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PCG::scalarSolve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
#6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#8  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#9  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
#10  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#11  ? at ??:?
#12  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#13  ? at ??:?
Thank for your help.

Last edited by Tibo99; May 28, 2020 at 18:13.
Tibo99 is offline   Reply With Quote

Old   May 28, 2020, 15:44
Default
  #4
Member
 
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2
superkelle is on a distinguished road
Maybe the problem is somewhere else, I mean the solver crashes for p calculation and need maximum iterations. And without the source everything works fine? Otherwise check BC.
superkelle is offline   Reply With Quote

Old   May 28, 2020, 18:22
Default
  #5
Member
 
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2
Tibo99 is on a distinguished road
Yes, everything work just fine if I don't apply the source term.

I see what you mean by the calculation for the pressure. The bounding value on epsilon are really high.....and the results on nut also, which is translate to a high value on k too(the picture is not there since I can't upload more than 5 files).

Since I use the following paper as a reference and a set of BC from my previous successful simulation without the source term, from what I understand, the BC don't need to be change because the source term on epsilon has been obtain by using the Inlet profile's equation shown in attachment and the BC are base from these profiles. The source term on k vanish. So there is no need to apply a source term on k with fvOptions.

Paper: https://www.researchgate.net/publica...n_of_ABL_flows
Attached Images
File Type: png formule.png (2.2 KB, 13 views)

Last edited by Tibo99; May 28, 2020 at 21:01.
Tibo99 is offline   Reply With Quote

Old   May 29, 2020, 01:44
Default
  #6
Member
 
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2
superkelle is on a distinguished road
Ok theres are now really just guesses...maybe first lower your relaxation factors? And you can try to lower yor source with a const factor like 1e-1 and look what happens. Do you have any restrictions on the mesh? Maybe refine it.
superkelle is offline   Reply With Quote

Old   May 29, 2020, 15:44
Default
  #7
Member
 
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2
Tibo99 is on a distinguished road
what I use for relaxation is: p=0.2, U=k=e=0.3, then I tried lowering these already and it didn't help. Still doing the same thing.

About multiplying by a factor, like I mentioned in the other thread that we talk together, I mutiplied the function by V[i], which is 1e-9 (cells volume=0.001mx0.001mx0.001m) and it did converge. So, it probably safe to say that the term V[i] has been use like a factor at this point.

What is limiting me on the mesh resolution is the resource, but I didn't try this option. The computer I use is not that powerful. For sure it is always better having a refined mesh but, Y+=50, which I think is relatively in the good range.

Thank again for your suggestion!

Stay safe!

Regards,

Last edited by Tibo99; May 30, 2020 at 14:20.
Tibo99 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Adding line source term for scalar transport by fvOptions wayne14 OpenFOAM Running, Solving & CFD 8 June 25, 2019 10:33
polynomial BC srv537 OpenFOAM Pre-Processing 4 December 3, 2016 09:07
[OpenFOAM.org] Error creating ParaView-4.1.0 OpenFOAM 2.3.0 tlcoons OpenFOAM Installation 13 April 20, 2016 17:34
SparceImage v1.7.x Issue on MAC OS X rcarmi OpenFOAM Installation 4 August 14, 2014 06:42
[swak4Foam] build problem swak4Foam OF 2.2.0 mcathela OpenFOAM Community Contributions 14 April 23, 2013 13:59


All times are GMT -4. The time now is 07:06.