
[Sponsors] 
May 20, 2019, 04:59 
From Laminar to Turbulent FLow

#1 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Hello Everyone,
I am using chtMultiRegionSimpleFoam and my Openfoam version is 4.1. I made my geometry in Salome and I am modelling a laminar fluid flow. Now I want to change the flow to Turbulent, But I don't know what changes do I need to make in my case to model a Turbulent flow. Can someone please guide me how to do this? Thank you 

May 20, 2019, 08:50 

#2 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
Hi. I am myself new to the platform but I can give you a few pointers.
Start with modifying the constant/turbulenceProperties where you set simulation type from laminar to RAS (or LES). Then there are a few keywords associated with the type of simulation you choose. eg. for RAS you have to specify the turbulence model used, eg kEpsilon.... you can easily search online for detailed list of available models and what are the corresponding keywords. Now, assuming you use kEpsilon model, the k and epsilon fields will be solved in order to obtain correct effective viscosity in our momentum equation (UEqn). So in your 0 directory, setup two files, named k and epsilon and define initial and boundary conditions. (You may have to declare k and omega if you use SSTkOmega model) These are all scalar quantities so you can start with the pressure file and modify it. Keep in mind to define dimensions etc properly. Alternatively, you can take a look at the high reynolds number cavity flow tutorial where you will find a sample of these files which you can change. Next, you will have to specify the following fluxes in your fvSchemes dictionary: div(phi,k), div(phi,epsilon) for which you can use Gauss Linear or Gauss vanLeer schemes as a starting point. This will tell the solver how these two scalars interact with the momentum flux in the discretized form of equations. Finally, you will need to setup numerical solution for the two discretized equations and error tolerance criteria. This is in the fvSolution dictionary where you may use the prescribed method for pressure and use it for the variables k and epsilon. You will be able to find information on this. Try to run this setup now. OF will tell you if it expects additional things to be declared, along with the dictionary where they are expected, which will give you an idea of where to make the changes. For example, you may need to define things like epsilonFinal which is the variable epsilon but used in the final iteration of a timestep (read about the PIMPLE iteration procedure to get an idea about this aspect of inner and outer iterations). Hope this helps. Note: at the time of running OF might ask you to also specify a third variable called nut which is not used by the simulation but seems necessary for the framework to work (I may be wrong here, this is just what I experienced in my case). 

May 20, 2019, 09:49 

#3 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Thank you so much for your reply.
I tell you what I have done. First, as you said, I changed my turbulenceProperties file from laminar to the following: turbulenceProperties Code:
simulationType RAS; RAS { RASModel kEpsilon; turbulence on; printCoeffs on; } And my changeDictionaryDict is given below: Code:
FoamFile { version 2.0; format ascii; class dictionary; object changeDictionaryDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // boundary { inlet { type patch; } outlet { type patch; } } U { internalField uniform (0 1e3 0); boundaryField { inlet { type fixedValue;//flowRateInletVelocity;//pressureInletVelocity;// //volumetricFlowRate 0.2; //extrapolateProfile yes; value uniform (0 2 0); } outlet { type inletOutlet;//zeroGradient;//// value $internalField;//uniform (0 0 0);// inletValue $internalField; } "fluid_to_box" { type noSlip; } } } T { internalField uniform 300; boundaryField { inlet { type fixedValue; value uniform 300;//$internalField; } outlet { type inletOutlet; value $internalField; inletValue $internalField; } "fluid_to_box" { type compressible::turbulentTemperatureCoupledBaffleMixed; Tnbr T; kappaMethod fluidThermo; value uniform 300; } } } epsilon { internalField uniform 0.01; boundaryField { inlet { type fixedValue; value uniform 0.01; } outlet { type zeroGradient;//inletOutlet; //inletValue uniform 0.01; } ".*" { type epsilonWallFunction; value uniform 0.01; } } } k { internalField uniform 0.1; boundaryField { inlet { type fixedValue;//inletOutlet; inletValue uniform 0.1; } outlet { type zeroGradient; //value uniform 0.1; } ".*" { type kqRWallFunction; value uniform 0.1; } } } p_rgh { internalField uniform 0; boundaryField { inlet { type fixedFluxPressure;//zeroGradient; value uniform 0; } outlet { type fixedValue; value uniform 0; } ".*" { type fixedFluxPressure; value uniform 0; } } } /*p { internalField uniform 0; boundaryField { ".*" { type zeroGradient; //value uniform 0; } } }*/ p { internalField uniform 0; boundaryField { inlet { type zeroGradient; //value uniform 1; } outlet { type fixedValue; value uniform 0; } "fluid_to_.*" { type zeroGradient; } } } Code:
ddtSchemes { default steadyState; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,U) bounded Gauss upwind; div(phi,K) bounded Gauss linear; div(phi,h) bounded Gauss upwind; div(phi,k) bounded Gauss linear;//bounded Gauss upwind; div(phi,K) bounded Gauss upwind; div(phi,epsilon) bounded Gauss linear;//bounded Gauss upwind; div(phi,R) bounded Gauss upwind; div(R) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default uncorrected; } fluxRequired { default no; p_rgh; } // ************************************************************************* // Code:
solvers { rho { solver PCG preconditioner DIC; tolerance 1e7; relTol 0; } p_rgh { solver GAMG; tolerance 1e7; relTol 0.01; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; maxIter 10; } "(Uhkepsilon)" { solver PBiCG; preconditioner DILU; tolerance 1e7; relTol 0.1; } } SIMPLE { momentumPredictor on; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 100000; rhoMin rhoMin [1 3 0 0 0] 0.2; rhoMax rhoMax [1 3 0 0 0] 2; } relaxationFactors { fields { rho 1; p_rgh 0.7; } equations { U 0.7; h 0.7; nuTilda 0.7; k 0.7; epsilon 0.7; omega 0.7; "ILambda.*" 0.7; } } // ************************************************************************* // Code:
Time = 0.1 Solving for fluid region fluid DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 0.06324153, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 0.01731178, No Iterations 2 DILUPBiCG: Solving for Uz, Initial residual = 1, Final residual = 0.05947683, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 0.9998487, Final residual = 0.01933637, No Iterations 2 Min/max T:300 300 GAMG: Solving for p_rgh, Initial residual = 1, Final residual = 0.7186809, No Iterations 10 time step continuity errors : sum local = 0.8787344, global = 0.03886439, cumulative = 0.03886439 Min/max rho:2 2 DILUPBiCG: Solving for epsilon, Initial residual = 0.2135582, Final residual = 0.01602102, No Iterations 1 bounding epsilon, min: 66.96574 max: 1520.873 average: 78.71354 DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 0.02653225, No Iterations 2 Solving for solid region box DICPCG: Solving for h, Initial residual = 0.8551727, Final residual = 0.03800197, No Iterations 2 Min/max T:300 300 Solving for solid region plate1 DICPCG: Solving for h, Initial residual = 1, Final residual = 0.009459444, No Iterations 2 Min/max T:300 300.0969 Solving for solid region plate2 DICPCG: Solving for h, Initial residual = 1, Final residual = 0.00875892, No Iterations 2 Min/max T:300 300.0972 Solving for solid region plate3 DICPCG: Solving for h, Initial residual = 1, Final residual = 0.008858153, No Iterations 2 Min/max T:300 300.1049 ExecutionTime = 1.56 s ClockTime = 1 s Time = 0.2 Solving for fluid region fluid DILUPBiCG: Solving for Ux, Initial residual = 0.6605907, Final residual = 0.02840858, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.6064321, Final residual = 0.0232098, No Iterations 1 DILUPBiCG: Solving for Uz, Initial residual = 0.5425838, Final residual = 0.0368132, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 0.1390125, Final residual = 0.003045702, No Iterations 2 Min/max T:3.050845e+07 3.376863e+07 GAMG: Solving for p_rgh, Initial residual = 0.1220398, Final residual = 0.0010412, No Iterations 7 time step continuity errors : sum local = 737166.4, global = 3599.479, cumulative = 3599.44 Min/max rho:2 2 DILUPBiCG: Solving for epsilon, Initial residual = 0.0009370397, Final residual = 3.433172e05, No Iterations 1 bounding epsilon, min: 3.263126e+09 max: 1.011055e+10 average: 3.31966e+07 DILUPBiCG: Solving for k, Initial residual = 0.426422, Final residual = 0.03064464, No Iterations 1 bounding k, min: 4.651673e+07 max: 1.451633e+08 average: 968806.9 Solving for solid region box DICPCG: Solving for h, Initial residual = 1, Final residual = 0.03428405, No Iterations 2 Min/max T:3173387 3195540 Solving for solid region plate1 DICPCG: Solving for h, Initial residual = 0.9985548, Final residual = 0.01910019, No Iterations 2 Min/max T:259.3801 682.9494 Solving for solid region plate2 DICPCG: Solving for h, Initial residual = 0.9995806, Final residual = 0.02288655, No Iterations 2 Min/max T:7493.021 3313.353 Solving for solid region plate3 DICPCG: Solving for h, Initial residual = 0.9994745, Final residual = 0.02052039, No Iterations 2 Min/max T:4433.842 12260.43 ExecutionTime = 2.37 s ClockTime = 2 s Time = 0.3 Solving for fluid region fluid DILUPBiCG: Solving for Ux, Initial residual = 0.9826942, Final residual = 0.0005316068, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.09963916, Final residual = 0.00187412, No Iterations 1 DILUPBiCG: Solving for Uz, Initial residual = 0.3857384, Final residual = 0.003870226, No Iterations 1 DILUPBiCG: Solving for h, Initial residual = 1, Final residual = 0.02362706, No Iterations 2 > FOAM FATAL ERROR: Maximum number of iterations exceeded From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::rhoConst<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy>] in file /home/ubuntu/OpenFOAM/OpenFOAM4.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66. FOAM aborting #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::error::abort() at ??:? #2 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::calculate() at ??:? #3 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::correct() at ??:? #4 ? at ??:? #5 __libc_start_main in "/lib/x86_64linuxgnu/libc.so.6" #6 ? at ??:? Aborted (core dumped) Thank you 

May 20, 2019, 14:39 

#4 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
I saw your query about your error on multiple posts. I am also curiously following the discussion.
In principle,the turbulence case setup correctly. However, in addition to the consistent error of max iterations exceeded, we can see that the value of k, epsilon are blowing up. In my case I had tackled this by setting boundary conditions correctly and it seems that the BCs you have set up are very loosely defined as almost all of them begin to get unbounded almost immediately. For starters: 1. Try to figure out the boundary conditions for k and epsilon at inlet by calculating approximate Reynolds number, turbulent mixing length, turbulence intesity ratio (There are prescriptions you can find for simple geometries which you can use.) 2. I am not too sure about the inletOutlet condition for velocity. Try a simpler definition like zeroGradient (atleast for the purpose of stability solution. You could try out better options later). In general you could replace a lot of BCs with simpler conditions like fixedValue of Gradient as these are primitive types of BCs in CFD. This is related to the advice you received about the pressurevelocity boundary conditions. It might also explain high value of residuals as system is unphysical. 3. I expect that with carefully configuring boundaries you will be able to make your solution stable, it takes some time but it is worthwhile to read about standard and suitable types for each variable. Now, coming to the final error, which seems to be related to species eqn, this post might be helpful : FOAM FATAL ERROR Maximum number of iterations exceeded ^but #3 can be expected to be heavily related to the other points as the divergence appears to propagate across variables before your solution crashes. These are only some of my hunches. Do let me know if you figure out what's happening. 

May 21, 2019, 03:12 

#5 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Thank you so much for your reply.
Here, I would like to add something. 1. Before I had one simple geometry(picture attached, narrow_pipe10), in which there is a rectangular fluid region, and the fluid flow was laminar. with the same boundary conditions. The solver was running upto steadystate. 2. Then, I just changed the shape of fluid region from rectangular to cylindrical, and put some curves in it (picture attached, inside_geometry). The three plates in both geometries are three heat sources which I made using fvOptions. My questions here are the following: 1. Does Laminar flow cannot be modeled in this curve shape fluid region? Should I change the flow from laminar to turbulent for this new geometry? because in this new geometry I am getting this error "maximum number of iterations exceeded", with the same boundary conditions as I had for the old geometry. 2. Can we model turbulent flow in chtMultiRegionSimpleFoam?, because I am using this solver. 3. How can I figure out the boundary conditions for k and epsilon using Reynolds number? (I am sorry for this question, as I am very new to Openfoam so I don't know much about it) These are my first basic confusion here, with same boundary condition one geometry is working while other is NOT. So, I thought may be I need to change the flow from Laminar to Turbulent for this new geometry. I shall be very thankful if you can clear these doubts. In the meanwhile, I look for the boundary conditions for k and epsilon using Reynolds number. Thank you 

May 21, 2019, 03:45 

#6 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
1. If the basic flow characteristics are similar (similar velocity, pressure drop per unit length etc), then geometry should not change things very drastically, except in regions of curvature etc. I cannot really say at this point whether it is due to a poor mesh, or boundary conditions need to be set better. Have you confirmed that the new geometry is meshed well? (It is anyway good practice to run checkMesh on your case)
You can make decision to switch to turbulence based on value of Reynolds number for your pipe crosssection. 2. From the file description it says it should be able to handle turbulence. (https://www.openfoam.com/documentati...8C_source.html) Code:
Description Steadystate solver for buoyant, turbulent fluid flow and solid heat conduction with conjugate heat transfer between solid and fluid regions. Code:
k = 3/2(UI)^2 I = I(Re) epsilon = C_mu( k^3/2)/l In your simplified heat exchanger, was your rectangular pipe also interacting with the heat sources? It seems to be located very differently for the purpose of heat exchange. 

May 21, 2019, 04:21 

#7 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Thank you for your prompt reply.
1. Yes, I have RUN the checkMesh utility to check my mesh. Here is the log for that: Code:
Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 31457 faces: 358864 internal faces: 355492 cells: 178589 faces per cell: 4 boundary patches: 3 point zones: 0 face zones: 0 cell zones: 5 Overall number of cells of each type: hexahedra: 0 prisms: 0 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 178589 polyhedra: 0 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces... Patch Faces Points Surface topology inlet 24 19 ok (nonclosed singly connected) outlet 24 19 ok (nonclosed singly connected) defaultFaces 3324 1674 ok (nonclosed singly connected) Checking geometry... Overall domain bounding box (0.02 0.07 0.059) (0.01 0.01 0.011) Mesh has 3 geometric (nonempty/wedge) directions (1 1 1) Mesh has 3 solution (nonempty) directions (1 1 1) Boundary openness (3.542138e18 1.54006e18 4.235165e19) OK. Max cell openness = 5.662842e16 OK. Max aspect ratio = 125.2255 OK. Minimum face area = 1.112013e08. Maximum face area = 1.179828e05. Face area magnitudes OK. Min volume = 1.387383e12. Max volume = 1.096243e08. Total volume = 4.2e05. Cell volumes OK. Mesh nonorthogonality Max: 88.07634 average: 21.58284 *Number of severely nonorthogonal (> 70 degrees) faces: 4086. Nonorthogonality check OK. <<Writing 4086 nonorthogonal faces to set nonOrthoFaces Face pyramids OK. Max skewness = 1.815118 OK. Coupled point location match (average 0) OK. Mesh OK. End Now, I am trying to calculate the Reynolds number of this new geometry. Here I have one question that, while calculating the Reynolds number, what would be the characteristic length? Is it the full length of pipe(from inlet to outlet) OR is it just the diameter of the pipe? Because if I take the full length of pipe then the Reynolds number turns out to be Turbulent, but If I select just the diameter of the pipe as a characteristic length, then it says the laminar. In my simplified geometry the pipe was not touching the heat sources, but there is a small gap in between. 

May 21, 2019, 04:55 

#8 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
There are quite a few elements with nonorthogonal faces, so that may be one of the areas to look into for potential improvement.
The diameter would be characteristic length so as to get an idea of flow through the crosssection. That is why I said changing geometry would not change things by much. Basically, I suppose you dont have to worry about turbulence at this point. The problem with your setup seems outside of this aspect. Well, one thing you could try is as follows: in changing the case you have changed the geometry of pipe and its heat exchange characteristics with the plates. Because you mention that your case setup works for the simple geometry (even though we may have scope to improve those BCs), it makes me think that if the problem comes because of temperature exchange at the fluidsolid contact and not because of the new mesh, then, if you used the simple rectangular pipe geometry, but brought it into contact with the walls, (similar to your actual case), then the simple geometry will also start to show similar error. This would mean that prime source of error is because the heat exchange interfaces are illdefined. Basically this way we can isolate the effect of new mesh, nonorthogonality, change in length etc and instead first perfect the simpler model which is working on some basic level, then it would be easier to take the setup to the other geometry. I would suggest to investigate that pipe geometry with similar contact as the final geometry, then you could perhaps see whether the error appears on change of pipe geometry, or due to the kind of heat exchange setup. If the original setup still runs to steady state, then it would mean that it is the new geometry which is causing divergence. In this case we may want to start looking at the nonorthogonal cells and try nNonOrthogonalCorrectors 1 to improve results. But move in small welldefined step changes that is very important. Do not make too many changes at once. It took me about a week or so to get my first simulation to become run, and that involved tweaking all these aspects. 

May 21, 2019, 05:57 

#9  
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Quote:
Hi... Following your suggestion, I made a new simple geometry, and it is running without errors, and as I increase the simulation time, it is converging to steadystate. I have attached picture of trials I have done at different simulation times (10s, 20s and 30s). you can see the temperature change in the fluid region due to heat sources, and After certain time, the whole fluid region and the plates take the same temperature. For your reference, I am attaching my boundary conditions here also: Code:
boundary { inlet { type patch; } outlet { type patch; } } U { internalField uniform (0 0 1e3); boundaryField { inlet { type fixedValue;//pressureInletVelocity;// //volumetricFlowRate 0.066; //extrapolateProfile yes; value uniform (0 0 0.2); } outlet { type zeroGradient; //value uniform (0 0 0);//$internalField; } /* ".*" { type noSlip; } */ "fluid_to_box" { type noSlip; } } } T { internalField uniform 300; boundaryField { inlet { type fixedValue; value $internalField; } outlet { type zeroGradient;//inletOutlet; value $internalField; //inletValue $internalField; } /*".*" { type zeroGradient; //value uniform 300; }*/ "fluid_to_box" { type compressible::turbulentTemperatureCoupledBaffleMixed; Tnbr T; kappaMethod fluidThermo; value uniform 300; } } } epsilon { internalField uniform 0.01; boundaryField { inlet { type fixedValue; value uniform 0.01; } outlet { type inletOutlet; inletValue uniform 0.01; } ".*" { type epsilonWallFunction; value uniform 0.01; } } } k { internalField uniform 0.1; boundaryField { inlet { type inletOutlet; inletValue uniform 0.1; } outlet { type zeroGradient; value uniform 0.1; } ".*" { type kqRWallFunction; value uniform 0.1; } } } p_rgh { internalField uniform 0; boundaryField { inlet { type zeroGradient; value uniform 0; } outlet { type fixedValue; value uniform 0; } ".*" { type fixedFluxPressure; value uniform 0; } } } p { internalField uniform 0; boundaryField { inlet { type zeroGradient;//fixedValue; //value uniform 10; } outlet { type fixedValue; value uniform 0; } "fluid_to_.*" { type zeroGradient; } } } // ************************************************************************* // 

May 21, 2019, 06:22 

#10 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
I am assuming this is a laminar solution as we have established that.
What I find strange is: is it normal to see the temperature of heat sources change as well? Anyway, so in this exact situation keeping everything else constant, if instead of the central rectangular pipe, we switch to a) cylindrical pipe (straight), and b) cylindrical pipe (S shaped curved) they all should technically run fine, as long as your direction of velocity flow is correctly defined (it should be moving into the inlet face  might be important to consider when using curved pipe). What do you think? 

May 21, 2019, 09:25 

#11  
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Quote:
Hi... As you suggested, I made one geometry with cylinder (image attached). and it is working fine. Then I changed the shape of cylindrical pipe. and then first it took a bit longer time in mesh generation. After that, running of simulation was also very slow, and then I got the following results.(image attached) Now, what do you think about these results? And one additional question here is that, how we check the direction of velocity?`By orientation of the geometry in Cartesian coordinate? for example, pipe is located in x direction, inlet is at (x=0) and outlet is at (x=10), then my velocity will be like this: type fixedValue value uniform (3 0 0); is it like this OR something else? 

May 21, 2019, 13:14 

#12  
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
Quote:


May 22, 2019, 03:12 

#13 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Hi...
I think the problem occurs when we try to model turbulent flow, because in curved shape pipe, I think it doesn't remain laminar anymore, and we modeled it as a laminar flow, that is why we are unable to see the expected results. OR do you think there could be some other problem? Now, the problem is how can we model turbulent flow in this curved shape pipe, and what are basics? Here is my understanding, please correct me if I am wrong: 1. We need to define a certain inlet velocity so that we can find the Reynolds number? 2. Once we have the Reynolds number, then we need to calculate the boundary conditions for k and epsilon? 3. And after that, we need to follow the steps that you mentioned in earlier posts relating to turbulent flows? kindly guide me, if I am missing something. Thank you 

May 22, 2019, 05:40 

#14 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
Why is the temperature in your curved pipe concentrated so asymmetrically? Have you double checked your direction of flow is correctly aligned in the new geometry? It is hard to tell based on your images. Can you check velocity vectors around inlet that it is oriented correctly?
If all that is correct then your curved pipe should also give a profile similar to the straight one if you have setup your inlet and outlet correctly and your new mesh is alright. I feel that if your velocity and pipe diameter are not changing much then neither should the Reynolds number, but you can/should confirm this about curved pipes, I do not know much about that aspect but it should be available in pipe and aerodynamics studies etc. 

May 22, 2019, 06:00 

#15  
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Quote:
Just an additional question here related to the flow direction, as you can see in the figure attached, the straight side of pipe is inlet and the curved side is outlet. And in the figure, you can see the Cartesian axis on the left bottom side. So, from this, can we say that our velocity at inlet is in the positive xdirection? that we can define like this: value uniform (velocity 0 0); 

May 22, 2019, 08:55 

#16  
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
Quote:
With this case, are you still getting the same error as your initial case, asked in first post? Otherwise it is quite strange if you actually got a solution with 1000 degrees at inlet and zero everywhere else, it seems very unphysical. In that case, the error had persisted with and without turbulence model, which is why I was concerned about the quality of BC and mesh setup. If your error is still there, it seems that either there is some flaw introduced due to the new mesh with curvature related issues in quality (remember 4000+ nonorthogonal faces?), or the physics of the case have changed due to tendencies of curvature in terms of flow and temperature. At this point, I am sorry that I can't precisely tell you what is happening, I am also curious to know. I am pretty sure it isn't a horrendous mistake we are making but a rather simple issue that needs to be fixed. See if playing with the mesh gets you somewhere. Also read about outerCorrectors and nonOrthogonalCorrectors they might be helpful in stabilizing your solution. Wherever in doubt, try replacing the BCs you don't fully understand with simple ones like fixedValue and zeroGradient. These initial divergences are typically mesh/BC problems. 

May 22, 2019, 10:45 

#17 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Thank you so much for your suggestion.
I will try to regenerate my mesh. In the meanwhile, I have one question related to Reynold's number. As we know the relation for Reynold's number: Re=(rho*u*L)/v where rho = density of the fluid (1000 in my case) u = velocity of fluid (1e3 in my case) L=Diameter of pipe(2.5mm in my case) v= kinematic Viscosity (959e6 in my case) Using the above values, my Reynolds number turns out to be less than 2300, it means that my flow will be Laminar. But if due to some reasons, the velocity of fluid changes somewhere in the pipe, would it change the fluid flow from Laminar to Turbulent? (because as per calculation using the above relation, if velocity increases more than 1.76m/s in my case, then the Re >4000, which means it is Turbulent flow) 

May 22, 2019, 11:59 

#18 
Member
Rishikesh
Join Date: Apr 2016
Posts: 63
Rep Power: 10 
Technically you are right, it will be a low turbulence situation. However, velocity going from 1e3 overall to 1.76 m/s in a local region isn't very common, that is a huge acceleration if you think about it  may occur if you change the crosssection suddenly by a large amount.
It may be possible if your average flow velocity itself is higher, then you'll have to use a turbulence model because your effective viscosity will have significant contribution of turbulence effect in addition to the normal laminar kinematic viscosity (959e6). This will change the flow patterns substantially and may even be necessary for convergence. 

May 22, 2019, 22:30 

#19 
Senior Member
Ruiyan Chen
Join Date: Jul 2016
Location: Hangzhou, China
Posts: 162
Rep Power: 9 
Hi Raza,
I think your problem is more related to the physics, not modelling (simulation). For example, is the fluid incompressible or compressible? If compressible, I assume that you are using some kind of equation of state to get density from temperature. I ask this question because in your output file, the maximum density reaches maxRho (2.0) really fast. Why did you bound your densities between 0.2 and 2.0? If the actual density of your fluid should be, say 5.0 (from some equation of state), but you are limiting it to be no larger than 2.0, of course your simulation blows up. Second, I myself is not familiar with chtMultiRegionSimpleFoam, but based on the description it is a steadystate solver. But, in your output file you are having Time = 0.1, 0.2, 0.3, which means you are not using deltaT = 1 in your controlDict. In this case of course your solver blows up, because you are using it to solve a transient case, which it is not designed for. Third, regarding the turbulence issues. The first thing you should do is to get a physical grasp of the type of the flow, mainly based on the Reynolds number. The curved geometry should not change the flow type, given your cross section area doesn't change. Curved geometry does give you a secondary flow, but that's not related to whether the flow being laminar or turbulent. With the geometry out of the way, here comes the real questionis your fluid incomressible or not? If it's compressible, its density changes with temperature, which means the density changes while heat transfer is happening. In this case, the density in the heated (or cooled) section might be very very different from what you have at the inlet. As a result, the flow changes from laminar to transition, and may eventually to turbulence. I don't really think this is what your case is about, but it's something to think about. Again, you should provide more information regarding the physics of your problem, the modelling or simulation part is secondary. I suggest the following: check if the fluid is incompressible or not. if it's incompressible, it would be much easier and we go from there; if it's compressible, what equation of state best describe your fluid, what about other thermophysical properties, such as viscosity? I feel like I asked a lot of questions, but hey, that's what you should do in CFD. Don't get lost in the modelling part of your problem (people do that a lot), start from fundamentals and physics. Figuring out how to condition your solver to match the physics is easy (although not as easy as you think), but figuring out the physics is much harder. Regards, Ruiyan 

May 23, 2019, 03:25 

#20 
Senior Member
Raza Javed
Join Date: Apr 2019
Location: Germany
Posts: 183
Rep Power: 7 
Hi Ruiyan,
Thank you so much for your reply. I have incompressible fluid in my simulation. And I have changed my controlDict file, you can see below: Code:
FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application chtMultiRegionSimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 100; deltaT 1; writeControl timeStep; writeInterval 1; purgeWrite 5; writeFormat ascii; writePrecision 7; writeCompression uncompressed; timeFormat general; timePrecision 6; runTimeModifiable true; // ************************************************************************* // But after this, my simulation stopped at T=35s. And the error is below: Code:
Time = 35 Solving for fluid region fluid DILUPBiCG: Solving for Ux, Initial residual = 0.6843924, Final residual = 0.02674027, No Iterations 2 DILUPBiCG: Solving for Uy, Initial residual = 0.6412737, Final residual = 0.03455567, No Iterations 2 DILUPBiCG: Solving for Uz, Initial residual = 0.5605676, Final residual = 0.02308348, No Iterations 2 DILUPBiCG: Solving for h, Initial residual = 0.6169205, Final residual = 0.01893477, No Iterations 3 > FOAM FATAL ERROR: Maximum number of iterations exceeded From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::rhoConst<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy>] in file /home/ubuntu/OpenFOAM/OpenFOAM4.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66. FOAM aborting #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::error::abort() at ??:? #2 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::calculate() at ??:? #3 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::correct() at ??:? #4 ? at ??:? #5 __libc_start_main in "/lib/x86_64linuxgnu/libc.so.6" #6 ? at ??:? Aborted (core dumped) I am also attaching the log file for your reference. 

Tags 
fluid, laminar, openfoam, turbulent 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Issues on the simulation of highspeed compressible flow within turbomachinery  dowlee  OpenFOAM Running, Solving & CFD  11  August 6, 2021 06:40 
Laminar and Turbulent flow of transient analysis  gbrajtm  Main CFD Forum  0  February 7, 2019 23:19 
Turbulent or laminar flow  yansheng  STARCCM+  4  July 1, 2016 17:53 
Do formula for Laminar and Turbulent Flow Calculation formulae change with Fluid  Parag Gadgil  FLUENT  0  June 19, 2012 07:07 
Can I use turbulent model to solve a laminar flow?  nikhil  FLUENT  5  February 1, 2011 10:42 