Help for the small implementation in turbulence model
Hi to all,
I am quite new about implementation and i want to change the kinematic eddy viscosity equation (nut) in k-epsilon turbulence model which is defined: Code:
nut_ = Cmu_*sqr(k_)/epsilon_; Code:
nut_ = f(rho)*Cmu_*sqr(k_)/epsilon_; Your help will be appreciated. Thanks in advance. B |
Hi,
maybe it'll be easier if you just multiply nut field returned from turbulence model by f(rho) inside the solver? |
Quote:
|
hi
Dear Alexey
Thanks so much for your advice. i have 2 questions. 1. When i look at the U.Eqn. in interphasechangeFoam: Code:
surfaceScalarField muEff 2. then i inserted my equation here as follows: Code:
surfaceScalarField muEff Code:
dency list for source file RNGinterPhaseChangeFoam.C thanks B |
1. mu is dynamic viscosity, nu is kinematic viscosity, mu = rho*nu. nut is kinematic turbulent viscosity.
2. Code:
fvc::interpolate(double Foam::pow((rho2-rho)/(rho2-rho1))*(rho1-rho2),5.0)+ rho2)*rho*turbulence->nut()) 2.2. You've got problems with brackets. According to your initial post the expressions should be Code:
fvc::interpolate((rho1 + (rho2 - rho)*Foam::pow(rho2 - rho1, -5))*rho*turbulence->nut()) Code:
(rho1 + (rho2 - rho)*Foam::pow(rho2 - rho1, -5)) Code:
(rho1.value() + (rho2.value() - rho.field())*Foam::pow((rho2 - rho1).value(), -5)) i.e. make your f(rho) non-dimensional. |
hi
Dear Alexey
Yeah you are right i just mixed with kinematic (nu) . and dynamic (mu) viscosities. As you know that rho1 and rho2 are constant. In order to avoid to mis write directly i wrote their values and defined the equation as follow( by the way at post #1 i wrote incorrect f(rho) true one is as follows): Code:
surfaceScalarField muEff Code:
Making dependency list for source file RNGinterPhaseChangeFoam.C Thanks in advance. B |
Unfortunately I don't have 2.1.1 installation, so these solutions will be quite theoretical (in 2.2.2 and 2.3.0 UEqn.H is slightly different).
As you need to make rho inside your f(rho) non-dimensional, you can introduce unit density densionedScalar and divide rho by this scalar: Code:
dimensionedScalar rhoUnit("rhoUnit", dimDensity, 1.0); Code:
- fvm::laplacian(muEff, U) Code:
volScalarField muEff("muEff", rho*nuEff()); |
Hi
Dear Alexey,
Thanks so much for your help. It is running now. Could you give some advice how can i improve myself about this kind of small implementations. On the other hand, T is transpose in the below equation?if not what is the meaning of it? fvc::div(muEff*dev(T(fvc::grad(U)))) Thanks in advance. B. |
Quote:
|
hi
Dear Alexey,
Thanks so much for your advices. I will have one another problem. I tried to use my same case in the another PC. (which has same version 2.1.1) When i run "wmake", the statement compile however it gives following warning: Code:
/opt/openfoam211/src/finiteVolume/lnInclude/readTimeControls.H: In function ‘int main(int, char**)’: HTML Code:
const bool adjustTimeStep = Could you tell me what is the mean reason of this problem? Thanks in advance. B |
Hi,
Cause you do not use directly maxDeltaT in your solver compiler complains about it. This warning appears during compilation of (almost) every OF solver. So you can ignore it. Concerning your question about decreasing time step, usually it means that you've got diverging solution. And there're plenty of reasons for it: inconsistent initial conditions, inconsistent boundary conditions, wrong settings in fvSolution, too large initial time step, etc. To tell you something about this problem I need to see the case files. |
hi
2 Attachment(s)
Dear Alexey,
First of all thank you so much for your instructive informations. Related with time decreasing problems i have attached the my initial conditions. First, i decreased deltaT then changed many times k and epsilon values but still it gives following error: HTML Code:
DILUPBiCG: Solving for alpha1, Initial residual = 0.00424273, Final residual = 1.84854e-06, No Iterations 1 thanks in advance. After checking the simulation i will also send some pictures. |
hi
4 Attachment(s)
Hello Alexey,
Also you can see the vel, liquid volume fraction and pres. distributions at attached pics. I really dont understand why it gives this kind of distribution. it is very strange that at the inlet pressure is increasing even mass flow (at the inlet) and outlet velocity seem constant???. And also, as you can see that cavitation takes place as attached the wall.. Could you tell me your idea. thanks in advance B. |
Hi,
1. Where is inlet, where is outlet? 2. Why do you think that 1 outer corrector step is enough for the simulation to converge? 3. Are you sure about 50 as IC for epsilon? |
Dear Alexey,
Thanks for reply. 1. Inlet is upstairs (long width) and down (small width) is outlet. 2. Actaully i dont know whether 1 outCorrecters enough or not? Could you advice what is the good setting? 2 or 3? 3. For epsilon i am not sure. According to initial calculation of epsilon based on k is 0.7 . I thought that it is very small. and before i read one post and one said that you can write 50~100 times higher for epsilon initial condition. So could you tell me your idea? Thanks in advance.. B |
Well,
1. As usual I'd suggest increasing nOuterCorrectors up to 50 and use residualControl to exit PIMPLE loop, i.e. change current PIMPLE dictionary in fvSolution from Code:
PIMPLE Code:
PIMPLE 2. Well, you can keep epsilon=50 if you like but if your simulation is blowing up why not try modifying this value? 3. I assume there's no problem with the mesh. Is it just simple rectangular mesh? Did you estimate y+? |
hi
Hi dear Alexey,
Thank so much for your reply. I started my calculation according to your advice. I will post results here when it finishes. I have some questions to make sure. 1. To increase the nOuterCorr means that we increase the calculation number for each time step right? and by adding "turbOnFinalIterOnly no" turbulent values are calculated in the same way for each time step? 2. yeah i changed the epsilon value and put calculated one as 0.7 this time. 3. yeah righ it is rectangular mesh and you can see the checkMesh results: HTML Code:
Mesh stats Thanks in advance. B |
If you take a look at the interPhaseChangeFoam source code, you'll find there:
Code:
while (pimple.loop()) - pimple.turbCorr() will return true on every iteration if turbOnFinalIterOnly is no, otherwise turbulence->correct() will be executed only once per time step. About y+: there are plenty of calculators. One of them in Tools menu - http://www.cfd-online.com/Tools/yplus.php. |
hi
Dear Alexey,
Finally, program gave again following floating point error:confused:. Mass flow, velocity inside of nozzle seem OK, but at the inlet pressure are increasing so excessively. Also inside of error i recognized that time step cont. error became crazy?? what is the meaning of this error and how can i correct it? Thanks in advance. B time step continuity errors : sum local = 4.64985e+125, global = -7.72928e+121, cumulative = -1.28369e+122 ERROR: HTML Code:
Courant Number mean: 1153.13 max: 617400 |
Hi,
well, now we arrived to the point where you can describe the problem you are trying to simulate ;) Cause 1e5 for p_rgh seems to be quite large. Also you can try to change relaxation settings: 0.7 for U, 0.3 for p_rgh, k and epsilon can be left as they are. |
Dear Alexey,
What do you mean exactly by saying "we arrived to the point where you can describe the problem you are trying to simulate" and how did you understand that 0.1 MPa is to large. In my simulation, outside is environment. Therefore, i set it atmospheric pressure. Ok, i will try 0 for pressure and will set the relaxation factors according to your advises. I will let you know results. Thanks so much B. |
Well,
I meant - "tell us about physical side of the problem you're trying to simulate", geometry, inflow/outflow conditions etc. Cause I have almost no knowledge about your problem, I can only guess what's going on and give general hints like "change relaxation factors, change PIMPLE dictionary" but the problem maybe elsewhere (for example in wrong direction of g). |
Hi
Dear Alexey,
I have tried following different cases: 1. 0/ k-epsilon I have changed wallFucntion setting =>as ZeroGradient 2. U inlet ==> FixedValue outlet= intletOutlet and outlet { type pressureInletOutletVelocity; value uniform (0 0 0); inletValue uniform (0 0 0); 3. i have used relaxation factors: Code:
relaxationFactors inlet ==> zeroGradient Outlet==> fixedValues ==>1e5, 1e4...0 (gradully decreased until 0) and also for outlet i tried: Code:
outlet All of these different challenges gives again following divergence error ::confused::confused::confused::confused: Code:
ime step continuity errors : sum local = 981070, global = -2.61762e-10, cumulative = 0.00804153 Thanks in advance. B |
Hi,
well, as you don't want to share an exact description of the physical problem you are trying to simulate, let's continue guessing. 1. You run modified version of interPhaseChangeFoam, have you tried running original solver? Just to try to localize what's causing troubles with convergence. 2. The only modification you've made was that scaling factor for nut? Or you've done something else? Can you share your modified code? I can continue with questions using just the pictures and solver output you've posted but think it won't be very productive. |
hi
Dear Alexey,
Thanks so much for your reply. I think that you misunderstood. Of course i can share exact description of the physical problem. I am trying to simulate cavitation inside of nozzle. My nozzle shape is as you can see at thread #13. 1. yeah i have tried with original interphasechangefoam. I had one thread before you can see here: http://www.cfd-online.com/Forums/ope...foam_help.html In original case, cavitation took place inside of nozzle, extended until some level and then disappeared. That was the problem. After that in some papers i saw that RANS turbulence viscosity is calculated high and need to have some modification to decrease the nut Now what i am trying to do. 2. Yes i have only modification which is bold in the file of U.Eqn as follow: Code:
dimensionedScalar rhoUnit("rhoUnit", dimDensity, 1.0); Please let me know what kind of info you want. Thanks B. |
Well, cause modeling of cavitation is not something that I do, my thoughts will be quite general (i.e. maybe completely useless ;)).
1. When you were doing simulations with original interPhaseChangeFoam, what was in your PIMPLE dictionary? Are you sure that the solution converged on every time step? 2. If you need lower turbulent viscosity why not start with just modifying Cmu (if you're using k-epsilon model). Set it 100 time lower (9e-4) and see what will happen. 3. Finally are you sure that settings for cavitation model are correct? In general you've got several points of possible errors: ICs/BCs, cavitation model settings, transport model settings, your modification of the code. See the influence of these things on the results, locate the source of the problem (well, I guess it'll be code modification as without it code was at least running and simulation wasn't blowing up). Plot muEff with your modification, compare it to the values without modifications. See if it even make any sense. |
All times are GMT -4. The time now is 08:50. |