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 Community New Posts Updated Threads Search

Like Tree4Likes

Reply
Page 1 of 2 1 2
 
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
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
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, 60 views)
File Type: jpg nut.jpg (101.1 KB, 33 views)
File Type: jpg Ux.jpg (96.5 KB, 37 views)
File Type: png Esource.png (3.3 KB, 59 views)
Attached Files
File Type: txt Epsilonterm.txt (26.3 KB, 22 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: 71
Rep Power: 6
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
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
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: 71
Rep Power: 6
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
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
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, 692 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: 71
Rep Power: 6
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
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
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

Old   June 6, 2020, 13:59
Default
  #8
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
Tibo99 is on a distinguished road
I revisited the finite volume method and I saw the following formula in attachment.

From what I saw in this formula, I was about right. The source term need to be multiply by the volume of the cell.
Attached Images
File Type: png fvm.png (24.5 KB, 71 views)
Tibo99 is offline   Reply With Quote

Old   June 6, 2020, 16:41
Default
  #9
HPE
Senior Member
 
HPE's Avatar
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 932
Rep Power: 12
HPE is on a distinguished road
Some of the new fvOptions related to the ABL modelling in the following link may help for the epsilon-/omega-based fvOptions you try to achieve? Have a look at `src/atmosphericModels/fvOptions/*` ?

https://develop.openfoam.com/Develop...e_requests/363

Hope this helps.
hwangpo and Tibo99 like this.
HPE is offline   Reply With Quote

Old   June 6, 2020, 18:08
Default
  #10
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
Tibo99 is on a distinguished road
Thank for the information!

I was trying to figure how the whole thing work with the 'fvOptions' framework by trying a simple case with a 'k' uniform profile.

Since I passed this step, my main goal has always been to implement an neutral ABL with a non-uniform 'k' and 'e' profile at the Inlet (varies with z). That give a bit more flexibility of matching wind tunnel results, which I want to achieve (for instance: https://www.sciencedirect.com/scienc...67610508001815).

That involve having new set of source term and BC.

Again, thank for the link. But, from what I've seen quickly is that the package is design for uniform 'k' profile which it don't apply in my case for the reason mentioned above. Even though, I always save good reference.

Best regards,
Tibo99 is offline   Reply With Quote

Old   June 6, 2020, 19:25
Default
  #11
HPE
Senior Member
 
HPE's Avatar
 
Herpes Free Engineer
Join Date: Sep 2019
Location: The Home Under The Ground with the Lost Boys
Posts: 932
Rep Power: 12
HPE is on a distinguished road
Hi,

Is there any reason why you think `k` is uniform? `k`-input can be spatiotemporal variant if it is a volScalarField or PatchFunction1<scalar> type. As far as I see, `k`-input of any of these functionalities is either of them. Do I miss something?

Thanks.
HPE is offline   Reply With Quote

Old   June 6, 2020, 21:21
Default
  #12
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
Tibo99 is on a distinguished road
From what I have noticed and understand in the literature is that using the k-e turbulence model for neutral ABL (steady-state), without any modification of the transport equation, the 'k' profile is uniform (spike can appear close to the wall though but become uniform after) and follow the relationship u*^2/Cmu^(1/2) and 'epsilon' u*^3/(kappa*(z+z0)). Specific BC at the Top and Bottom must be added if you want to match the velocity profile U obtain at the Outlet compare to the one used at the Inlet. (see: https://onlinelibrary.wiley.com/doi/....1002/fld.2709)

If I want to change the 'k' profile from uniform to a custom profile/log profile at the Inlet, often used to match wind tunnel results (see the paper I linked in the previous post I did), to obtain a neutral ABL and having the corresponding results at the Outlet, there is 2 options in order to reflecting these new conditions at the Inlet:

1- Changing the transport equation,
2- Adding a source term to the original transport equation (see: https://www.researchgate.net/publica...n_of_ABL_flows).

Obviously, the BC at the Top and Bottom must be different also.

So, what I have seen from the results presented for 'k' in the link you shared is that you can see what we call the 'spike' and then it become basically uniform. So, this make me think that the transport equation, the fvOptions and the wall functions are not adapted for a new type of 'k' profile that I'm looking for at the Inlet.
HPE likes this.
Tibo99 is offline   Reply With Quote

Old   June 7, 2020, 11:14
Default
  #13
Senior Member
 
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6
Tibo99 is on a distinguished road
Hi @Nasim99!

No problem. If you find something new or if you find something that I said is wrong and you can confirm it, don't hesitate to comment.

If you try to do something similar, it would be nice to compare and share the BC. It don't converge right now and since I confirmed that I use the right source term for my case, I think the issue is from the BC.

The only thing I wondering is, if it's possible to use Cmu constant if 'k' is non-uniform(varies with 'z') and adding a source term instead of modifying the transport equation? If someone could confirm that, it would be appreciated.

Regards,
Last edited by Tibo99; June 7, 2020 at 23:58.
Tibo99 is offline   Reply With Quote

Old   June 30, 2021, 12:22
Default
  #14
Member
 
Kabir Shariff
Join Date: Oct 2016
Location: France
Posts: 53
Rep Power: 9
Kbshariff is on a distinguished road
Hello,

I want to apply a source term to the k-epsilon model using the fvOptions ScalarcodedSource. Attached herewith is the source term I want to apply from the literature.

This source term is applied in a particular location defined as cellZone in topoSet
(see attached schematic)

I first applied the source term to the k field only to verify the code is ok before adding the epsilon source.

my problem is that
The k & epsilon becomes zero just after few cells from the inlet (see attached file kSource.csv) but the fields gives good result without the source terms

the source term is applied on the cellzone zone1 as shown in the logfile below.
how can I rectify the problem please? is my code ok?

attached herewith is my fvOptions code for the kSource

Code:
kSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone         zone1;
    fields          (k);
    name            kSource;

        codeCorrect
        #{
        #};

    codeConstrain
    #{
    #};

    codeAddSup
    #{
				
        const scalar epsIn_ = 1.38e-4;  // 5%

        const vectorField& CellC = mesh_.C(); 
        
		
	scalarField& kSource = eqn.source();
				
        // Apply the source

        // Only apply when we have reached the start time

        forAll(CellC, i)
        {
          							
               kSource[i] += epsIn_;
            			
        };
        
        
    #};
}
the log file
Code:
Selecting finite volume options type actuationDiskSource
    Source: disk1
    - selecting cells using cellSet actuationDisk1
    - selected 1552 cell(s) with volume 0.00596923076923
    - selecting cells using cellSet actuationDisk1
    - creating actuation disk zone: disk1
    - force computation method: Froude
Selecting finite volume options type scalarCodedSource
    Source: kSource
    - selecting cells using cellZone zone1
    - selected 32980 cell(s) with volume 0.198724573946

Starting time loop

turbulenceFields turbulenceFields: storing fields:
    turbulenceProperties:I

    weight field = none

surfaceFieldValue Uupstream:
    operation     = weightedAreaAverage
    weight field  = none


Time = 1

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0449048902513, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.999982619441, Final residual = 0.024345582838, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 1, Final residual = 0.0884208859647, No Iterations 2
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.000968185944403, No Iterations 124
time step continuity errors : sum local = 4.3030486421e-05, global = -2.0809324322e-05, cumulative = -2.0809324322e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.149887418611, Final residual = 0.0111000171302, No Iterations 2
Using dynamicCode for fvOption::kSource at line 37 in "/dlocal/run/8723855/constant/fvOptions.kSource"
Could not load "/dlocal/run/8723855/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libkSource_9ec1fe8a1a78c2055ca6aaec53cc716f66c67a36.so"
/dlocal/run/8723855/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libkSource_9ec1fe8a1a78c2055ca6aaec53cc716f66c67a36.so: cannot open shared object file: No such file or directory
Creating new library in "dynamicCode/kSource/platforms/linux64GccDPInt32Opt/lib/libkSource_9ec1fe8a1a78c2055ca6aaec53cc716f66c67a36.so"
Invoking wmake libso /dlocal/run/8723855/dynamicCode/kSource
wmake libso /dlocal/run/8723855/dynamicCode/kSource
    ln: ./lnInclude
    dep: codedFvOptionTemplate.C
    Ctoo: codedFvOptionTemplate.C
    ld: /dlocal/run/8723855/dynamicCode/kSource/../platforms/linux64GccDPInt32Opt/lib/libkSource_9ec1fe8a1a78c2055ca6aaec53cc716f66c67a36.so
Selecting finite volume options type kSource
    Source: kSource
    - selecting cells using cellZone zone1
    - selected 32980 cell(s) with volume 0.198724573946
smoothSolver:  Solving for k, Initial residual = 0.999999999999, Final residual = 0.0825490627951, No Iterations 3
bounding k, min: -0.350071661333 max: 0.0575828573434 average: -0.0708006843692
ExecutionTime = 15.84 s  ClockTime = 20 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.871096205218, Final residual = 0.0145172664649, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.4664537379, Final residual = 0.0126773675555, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.495367135106, Final residual = 0.042222418612, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.0194949736121, Final residual = 1.7095166458e-05, No Iterations 108
time step continuity errors : sum local = 0.000309894273643, global = 9.41105615543e-05, cumulative = 7.33012372323e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.59019740074, Final residual = 4.25091897642e-13, No Iterations 1
bounding epsilon, min: 2.84968168995e-21 max: 0.0383230234325 average: 0.000116078627665
smoothSolver:  Solving for k, Initial residual = 7.39851510179e-08, Final residual = 7.39851510179e-08, No Iterations 0
ExecutionTime = 28.72 s  ClockTime = 33 s

Time = 3

smoothSolver:  Solving for Ux, Initial residual = 0.367516183527, Final residual = 0.00927443461916, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.480053536097, Final residual = 0.0305014083565, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 0.6118937401, Final residual = 0.0493583388814, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.0600169917322, Final residual = 5.30007958956e-05, No Iterations 84
time step continuity errors : sum local = 0.000283702185864, global = 9.06862878832e-05, cumulative = 0.000163987525116
smoothSolver:  Solving for epsilon, Initial residual = 0.565279741295, Final residual = 3.35983386159e-11, No Iterations 1
bounding epsilon, min: 2.84968168995e-21 max: 0.502567369362 average: 0.000383204117346
smoothSolver:  Solving for k, Initial residual = 2.49344628981e-07, Final residual = 2.49344628981e-07, No Iterations 0
ExecutionTime = 39.07 s  ClockTime = 43 s
Attached Images
File Type: png el-kasmi.PNG (8.6 KB, 11 views)
File Type: jpg domain-elkasmi.JPG (22.2 KB, 11 views)
Attached Files
File Type: txt ksource.txt (34.5 KB, 7 views)
Kbshariff is offline   Reply With Quote

Old   July 1, 2021, 08:29
Default
  #15
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
I think you have to use the correct cell indices like it is done here:

https://github.com/OpenFOAM/OpenFOAM...emplates.C#L62

Code:
kSource[CellC[i]] += epsIn_;
jherb is offline   Reply With Quote

Old   July 1, 2021, 11:27
Default
  #16
Member
 
Kabir Shariff
Join Date: Oct 2016
Location: France
Posts: 53
Rep Power: 9
Kbshariff is on a distinguished road
Quote:
Originally Posted by jherb View Post
I think you have to use the correct cell indices like it is done here:

https://github.com/OpenFOAM/OpenFOAM...emplates.C#L62

Code:
kSource[CellC[i]] += epsIn_;
Thank you,

I made the correction as suggested but I get this error message

fvOptions
Code:
kSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone         zone1;
    fields          (k);
    name            kSource;

        codeCorrect
        #{
        #};

    codeConstrain
    #{
    #};

    codeAddSup
    #{
				
        const scalar epsIn_ = 1.38e-4;  // 5%

        const vectorField& CellC = mesh_.C(); 
        
		
	scalarField& kSource = eqn.source();
				
        // Apply the source

        // Only apply when we have reached the start time

        forAll(CellC, i)
        {
          							
               kSource[CellC[i]] += epsIn_;
            			
        };
        
        
    #};
}
logfile
Code:
Creating finite volume options from "constant/fvOptions"

Selecting finite volume options type scalarCodedSource
    Source: kSource
    - selecting cells using cellZone zone1
    - selected 32980 cell(s) with volume 0.198724573946

Starting time loop

turbulenceFields turbulenceFields: storing fields:
    turbulenceProperties:I

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0449048902513, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.999982619441, Final residual = 0.024345582838, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 1, Final residual = 0.0884208859647, No Iterations 2
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.000968185944403, No Iterations 124
time step continuity errors : sum local = 4.3030486421e-05, global = -2.0809324322e-05, cumulative = -2.0809324322e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.149887418611, Final residual = 0.0111000171302, No Iterations 2
Using dynamicCode for fvOption::kSource at line 21 in "/dlocal/run/8724992/constant/fvOptions.kSource"
Could not load "/dlocal/run/8724992/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libkSource_57ba7fa6589edcb86f349c030adb191d5977367b.so"
/dlocal/run/8724992/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libkSource_57ba7fa6589edcb86f349c030adb191d5977367b.so: cannot open shared object file: No such file or directory
Creating new library in "dynamicCode/kSource/platforms/linux64GccDPInt32Opt/lib/libkSource_57ba7fa6589edcb86f349c030adb191d5977367b.so"
Invoking wmake libso /dlocal/run/8724992/dynamicCode/kSource
wmake libso /dlocal/run/8724992/dynamicCode/kSource
    ln: ./lnInclude
    dep: codedFvOptionTemplate.C
    Ctoo: codedFvOptionTemplate.C
/dlocal/run/8724992/constant/fvOptions.kSource: In member function ‘virtual void Foam::fv::kSourceFvOptionscalarSource::addSup(Foam::fvMatrix<double>&, Foam::label)’:
/dlocal/run/8724992/constant/fvOptions.kSource:51:23: error: no match for ‘operator[]’ (operand types are ‘Foam::scalarField {aka Foam::Field<double>}’ and ‘const Foam::Vector<double>’)
/dlocal/run/8724992/constant/fvOptions.kSource:51:23: note: candidates are:
In file included from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/UList.H:652:0,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/List.H:46,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/HashTable.C:33,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/Istream.H:258,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/ISstream.H:42,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/IOstreams.H:40,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/VectorSpace.C:29,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/VectorSpace.H:279,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/Vector.H:48,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/vector.H:42,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/fieldTypes.H:37,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvMatricesFwd.H:34,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/finiteVolume/lnInclude/fvOption.H:50,
                 from /soft/OpenFOAM-v2006/OpenFOAM-v2006/src/fvOptions/lnInclude/cellSetOption.H:57,
                 from codedFvOptionTemplate.H:105,
                 from codedFvOptionTemplate.C:29:
/soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/UListI.H:246:11: note: T& Foam::UList<T>::operator[](Foam::label) [with T = double; Foam::label = int]
 inline T& Foam::UList<T>::operator[](const label i)
           ^
/soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/UListI.H:246:11: note:   no known conversion for argument 1 from ‘const Foam::Vector<double>’ to ‘Foam::label {aka int}’
/soft/OpenFOAM-v2006/OpenFOAM-v2006/src/OpenFOAM/lnInclude/UListI.H:256:17: note: const T& Foam::UList<T>::operator[](Foam::label) const [with T = double; Foam::label = int]
 inline const T& Foam::UList<T>::operator[](const label i) const
                 ^
Kbshariff is offline   Reply With Quote

Old   July 1, 2021, 12:12
Default
  #17
Senior Member
 
Joachim Herb
Join Date: Sep 2010
Posts: 650
Rep Power: 21
jherb is on a distinguished road
I checked again: https://www.openfoam.com/documentati...ces-coded.html


mesh_.C() seems to be cell centres cooridnates. You want the list of the cells in the cell set. See the example in this bug report: https://bugs.openfoam.org/view.php?id=2361

Modified with your source:
Code:
                const labelList& cellIDs = cells();

                forAll(cellIDs, i)
                {
                    label cellI = cellIDs[i];
                    kSource[cellI] += epsIn_;
                }
How did I find it: Googled for coded source selectionmode cellzone openfoam
jherb is offline   Reply With Quote

Old   July 2, 2021, 12:37
Default
  #18
Member
 
Kabir Shariff
Join Date: Oct 2016
Location: France
Posts: 53
Rep Power: 9
Kbshariff is on a distinguished road
Quote:
Originally Posted by jherb View Post
I checked again: https://www.openfoam.com/documentati...ces-coded.html


mesh_.C() seems to be cell centres cooridnates. You want the list of the cells in the cell set. See the example in this bug report: https://bugs.openfoam.org/view.php?id=2361

Modified with your source:
Code:
                const labelList& cellIDs = cells();

                forAll(cellIDs, i)
                {
                    label cellI = cellIDs[i];
                    kSource[cellI] += epsIn_;
                }
How did I find it: Googled for coded source selectionmode cellzone openfoam

Thank you very much for the assistance.

the simulation is running but the result is sketchy


Code:
kSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone         zone1;
    fields          (k);
    name            kSource;

        codeCorrect
        #{
        #};

    codeConstrain
    #{
    #};

    codeAddSup
            #{

                scalarField& keSource = eqn.source();
       		const scalar epsIn_ = 1.38e-4;
                const labelList& cellIDs = cells();

                forAll(cellIDs, i)
                {
                    label cellI = cellIDs[i];
                    keSource[cellI] -= epsIn_;
                }
            #};
}

epsilonSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone         zone1;
    fields          (epsilon);
    name            epsilonSource;

        codeCorrect
        #{
        #};

    codeConstrain
    #{
    #};

    codeAddSup
            #{

                scalarField& epsilonSource = eqn.source();
       		const scalar epsIn_	= 1.38e-4;
		const scalar kIn_ 	= 0.0024;
		const scalar C_		= 1.92;
                const labelList& cellIDs = cells();

                forAll(cellIDs, i)
                {
                    label cellI = cellIDs[i];
                    epsilonSource[cellI] -= C_*epsIn_/kIn_;
                }
            #};
}
the source term is multiplied by 'rho' in the literature. I did not use it since rho is unity by default.

The results show the source term added in k & epsilon in the region but reduces to zero afterwards. I attached herewith a comparison between source and no source.

Thank you very much. Have a nice weekend
Attached Images
File Type: jpg TI.jpg (19.0 KB, 17 views)
File Type: jpg k.jpg (19.0 KB, 10 views)
File Type: jpg epsilon.jpg (16.8 KB, 13 views)
Kbshariff is offline   Reply With Quote

Old   February 9, 2023, 12:22
Default
  #19
Lil
New Member
 
Lila
Join Date: Nov 2019
Posts: 7
Rep Power: 6
Lil is on a distinguished road
Hello,


Can you please tell me what is "epsilonSource()" in the kepsilon.C file ?
Lil is offline   Reply With Quote

Old   February 10, 2023, 06:28
Default
  #20
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 670
Rep Power: 14
Tobermory will become famous soon enough
Take a look at kEpsilon.H and you will see that it is a virtual function for the source term in the epsilon equation:
Code:
         virtual tmp<fvScalarMatrix> epsilonSource() const;
EDIT: if you are still struggling, take a look at the buoyantKEpsilon model, which is derived from the kEpsilon class, and you will see the kSource and epsilonSource functions in action.
Last edited by Tobermory; February 12, 2023 at 12:27.
Tobermory is offline   Reply With Quote

Reply
Page 1 of 2 1 2

« Previous Thread | Next Thread »

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 17:25.

Contact Us - CFD Online - Privacy Statement - Top