
[Sponsors] 
Adding source term with fvOptions to ke model 

LinkBack  Thread Tools  Search this Thread  Display Modes 
May 21, 2020, 16:34 
Adding source term with fvOptions to ke model

#1 
Member
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2 
Hi everyone!
I got some issue to implement a source term via fvOptions to the ke model. The source term is apply on epsilon(only) to the ke 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, ycoordinate 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 epsilonBC 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.3474e05; 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 #{ #}; } // ************************************************************************* // 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 : reneOptiPlex790 PID : 16903 I/O : uncollated Case : /home/rene/OpenFOAM/renev1906/run/1D_E1 nProcs : 1 trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE). fileModificationChecking : Monitoring runtime modified files using timeStampMaster (fileModificationSkew 10) allowSystemOperations : Allowing usersupplied system call operations // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time Create mesh for time = 0 SIMPLE: convergence criteria field U tolerance 1e06 field k tolerance 1e06 field epsilon tolerance 1e06 field omega tolerance 1e06 field nuTilda tolerance 1e06 field T tolerance 1e06 field p_rgh tolerance 1e06 field p tolerance 1e06 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 highRe 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.5e07 Starting time loop turbulenceFields turbulenceFields1: storing fields: turbulenceProperties:R turbulenceProperties:I turbulenceProperties:L Time = 1 smoothSolver: Solving for Ux, Initial residual = 5.781665123e06, Final residual = 2.320453622e10, No Iterations 8 smoothSolver: Solving for Uy, Initial residual = 0.9089618075, Final residual = 5.582492529e05, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 1.786751413e15, Final residual = 4.823728453e16, No Iterations 1 GAMG: Solving for p, Initial residual = 0.7530448709, Final residual = 5.369605972e05, No Iterations 21 time step continuity errors : sum local = 3.673528501e19, global = 3.673528501e19, cumulative = 3.673528501e19 Using dynamicCode for fvOption::epsilonSource at line 1 in "/home/rene/OpenFOAM/renev1906/run/1D_E1/system/fvOptions.epsilonSource" Creating new library in "dynamicCode/epsilonSource/platforms/linux64GccDPInt32Opt/lib/libepsilonSource_273d0b6b78f639d7c93a407718428cb1cf5e4e35.so" Invoking wmake libso /home/rene/OpenFOAM/renev1906/run/1D_E1/dynamicCode/epsilonSource wmake libso /home/rene/OpenFOAM/renev1906/run/1D_E1/dynamicCode/epsilonSource dep: codedFvOptionTemplate.C Ctoo: codedFvOptionTemplate.C ld: /home/rene/OpenFOAM/renev1906/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.5e07 smoothSolver: Solving for epsilon, Initial residual = 0.9998814658, Final residual = 2.827019512e05, No Iterations 8 bounding epsilon, min: 506399185.4 max: 2140.790268 average: 5800585.107 smoothSolver: Solving for k, Initial residual = 1, Final residual = 9.257761502e05, No Iterations 7 ExecutionTime = 0.06 s ClockTime = 3 s Time = 2 smoothSolver: Solving for Ux, Initial residual = 0.0002397693829, Final residual = 8.113253956e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.01029908582, Final residual = 4.358351871e07, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.0629515603, Final residual = 2.569515862e06, 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.213998486e05, global = 1.400000834e09, cumulative = 1.400000834e09 smoothSolver: Solving for epsilon, Initial residual = 0.003063504061, Final residual = 2.350539837e07, No Iterations 6 bounding epsilon, min: 4.792808503e11 max: 1947.869938 average: 13.85865537 smoothSolver: Solving for k, Initial residual = 0.397705588, Final residual = 2.981436582e05, No Iterations 7 ExecutionTime = 0.3 s ClockTime = 3 s Time = 3 smoothSolver: Solving for Ux, Initial residual = 6.542876533e05, Final residual = 3.414747811e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.008345991815, Final residual = 6.916702663e07, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.4207610645, Final residual = 1.767895306e05, No Iterations 7 GAMG: Solving for p, Initial residual = 0.9875079438, Final residual = 2.876784968e05, No Iterations 5 time step continuity errors : sum local = 72.2530223, global = 7.186872358e10, cumulative = 2.11868807e09 smoothSolver: Solving for epsilon, Initial residual = 1.779790796e10, Final residual = 2.332839125e11, 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.109144553e05, No Iterations 7 ExecutionTime = 0.3 s ClockTime = 3 s Time = 4 smoothSolver: Solving for Ux, Initial residual = 0.0001047864544, Final residual = 3.802370649e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 1.067786309e09, Final residual = 9.803739005e11, No Iterations 1 smoothSolver: Solving for Uz, Initial residual = 5.08984362e07, Final residual = 3.922613063e11, No Iterations 6 GAMG: Solving for p, Initial residual = 0.9999368938, Final residual = 9.488941883e05, No Iterations 974 time step continuity errors : sum local = 94.93563983, global = 2.048071334e05, cumulative = 2.048283203e05 smoothSolver: Solving for epsilon, Initial residual = 1.464161316e10, Final residual = 2.944491158e11, No Iterations 1 bounding epsilon, min: 2.276798111e11 max: 1591232665 average: 14897422.52 smoothSolver: Solving for k, Initial residual = 4.594718889e09, Final residual = 8.271137915e11, No Iterations 3 ExecutionTime = 0.36 s ClockTime = 3 s Time = 5 smoothSolver: Solving for Ux, Initial residual = 5.71021682e05, Final residual = 3.086792287e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 1.395646902e09, Final residual = 6.012380828e11, No Iterations 2 smoothSolver: Solving for Uz, Initial residual = 1.899201834e10, Final residual = 3.854897013e11, No Iterations 1 GAMG: Solving for p, Initial residual = 0.9998158918, Final residual = 6.610612207e05, No Iterations 10 time step continuity errors : sum local = 15182.47696, global = 0.005463523176, cumulative = 0.005443040344 smoothSolver: Solving for epsilon, Initial residual = 5.515123834e09, Final residual = 3.612469016e11, No Iterations 4 bounding epsilon, min: 1.555379741e11 max: 3.206587777e+12 average: 2.633791325e+10 smoothSolver: Solving for k, Initial residual = 8.865456493e07, Final residual = 6.354148454e11, No Iterations 7 ExecutionTime = 0.37 s ClockTime = 3 s Time = 6 smoothSolver: Solving for Ux, Initial residual = 4.903075674e05, Final residual = 2.76209214e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.001152352045, Final residual = 3.81812229e08, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.00082024583, Final residual = 5.758955053e08, No Iterations 6 GAMG: Solving for p, Initial residual = 0.9912377818, Final residual = 6.390959052e05, 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.672863807e08, No Iterations 7 bounding epsilon, min: 1.131701388e11 max: 7.930001666e+12 average: 8.597395232e+10 smoothSolver: Solving for k, Initial residual = 0.01084036878, Final residual = 7.70259312e07, No Iterations 7 ExecutionTime = 0.38 s ClockTime = 3 s Time = 7 smoothSolver: Solving for Ux, Initial residual = 0.0001017155284, Final residual = 4.623390674e09, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.005494247008, Final residual = 1.851745399e07, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.01430007137, Final residual = 4.569324173e07, No Iterations 6 GAMG: Solving for p, Initial residual = 0.8459483358, Final residual = 3.237544819e05, No Iterations 12 time step continuity errors : sum local = 6614.508604, global = 13.01227121, cumulative = 196.0737873 smoothSolver: Solving for epsilon, Initial residual = 7.049825885e06, Final residual = 3.746246606e10, No Iterations 7 bounding epsilon, min: 7.098149768e12 max: 1.287364284e+14 average: 9.485271478e+11 smoothSolver: Solving for k, Initial residual = 0.001794322734, Final residual = 1.348391127e07, No Iterations 7 ExecutionTime = 0.38 s ClockTime = 3 s Time = 8 smoothSolver: Solving for Ux, Initial residual = 0.001637598567, Final residual = 6.434227911e08, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.1603902233, Final residual = 4.344377787e06, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 0.0002409532143, Final residual = 1.363044774e08, 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.42542935e11, Final residual = 2.006210096e11, No Iterations 1 smoothSolver: Solving for k, Initial residual = 0.4804684285, Final residual = 4.048431869e05, No Iterations 7 ExecutionTime = 0.59 s ClockTime = 4 s Time = 9 smoothSolver: Solving for Ux, Initial residual = 0.000961475335, Final residual = 4.019731076e08, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.05898306748, Final residual = 1.477984489e06, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 3.178084763e08, Final residual = 2.804660479e11, 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.413948522e17, Final residual = 7.529783842e18, No Iterations 1 smoothSolver: Solving for k, Initial residual = 0.2539284515, Final residual = 2.201685003e05, No Iterations 7 ExecutionTime = 0.84 s ClockTime = 4 s Time = 10 smoothSolver: Solving for Ux, Initial residual = 0.0009077776947, Final residual = 3.482148807e08, No Iterations 7 smoothSolver: Solving for Uy, Initial residual = 0.1731602652, Final residual = 7.857890645e06, No Iterations 7 smoothSolver: Solving for Uz, Initial residual = 2.73445177e08, Final residual = 2.887628211e11, 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.155609824e11, No Iterations 1 smoothSolver: Solving for k, Initial residual = 1, Final residual = 9.253973429e05, No Iterations 7 ExecutionTime = 1.33 s ClockTime = 4 s Time = 11 smoothSolver: Solving for Ux, Initial residual = 1.40996752e13, Final residual = 5.783231539e15, No Iterations 1 smoothSolver: Solving for Uy, Initial residual = 9.056049294e07, Final residual = 4.438814615e11, No Iterations 6 smoothSolver: Solving for Uz, Initial residual = 1.40996752e13, Final residual = 5.783231539e15, No Iterations 1 Last edited by Tibo99; May 24, 2020 at 10:22. 

May 28, 2020, 07:09 

#2 
Member
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2 
Hi, do you get an floating point error after t=11? I am not sure about the "+=" maybe try simple "=".


May 28, 2020, 12:49 

#3 
Member
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2 
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.40996752e13, Final residual = 5.783231539e15, No Iterations 1 smoothSolver: Solving for Uy, Initial residual = 9.056049294e07, Final residual = 4.438814615e11, No Iterations 6 smoothSolver: Solving for Uz, Initial residual = 1.40996752e13, Final residual = 5.783231539e15, No Iterations 1 #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 ? in /lib/x86_64linuxgnu/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_64linuxgnu/libc.so.6 #13 ? at ??:? Last edited by Tibo99; May 28, 2020 at 18:13. 

May 28, 2020, 15:44 

#4 
Member
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2 
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.


May 28, 2020, 18:22 

#5 
Member
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2 
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 Last edited by Tibo99; May 28, 2020 at 21:01. 

May 29, 2020, 01:44 

#6 
Member
alexander thierfelder
Join Date: Dec 2019
Posts: 41
Rep Power: 2 
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 1e1 and look what happens. Do you have any restrictions on the mesh? Maybe refine it.


May 29, 2020, 15:44 

#7 
Member
René Thibault
Join Date: Dec 2019
Posts: 32
Rep Power: 2 
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 1e9 (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. 

Thread Tools  Search this Thread 
Display Modes  


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 PreProcessing  4  December 3, 2016 09:07 
[OpenFOAM.org] Error creating ParaView4.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 