Implementing new turbulence model
2 Attachment(s)
I want to implement a new turbulence model (DafaliasYounis) into OpenFOAM. This turbulence model is a Reynolds Stress model, which models the pressurestrain correlation Phi as shown in equation (10) attached to this post. Sij, Wij and bij are defined as in the second attached jpg.
I am working on this implementation now for quite some time and I encountered very much problems. Fortunately hjasak, a member of this forum, helped me a lot (Thanks again). Thus, I will post here the messages we exchanged, so that it is possible to get a full overview about the problems we discussed. I still could not implement the model, so I am very greatful for every help. Thanks a lot! 
Dear Mr Jasak,
I am a student of Aerospace Engineering, doing my master thesis at the University of California in Davis. For my thesis, I try to implement a new turbulence model in OpenFoam. For doing that, I want to edit the already existing LaunderGibsonRSTM. I edited the original equation (LaunderGibsonRSTM.C file) like this: Code:
tmp<fvSymmTensorMatrix> REqn Code:
+C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_)))&b) http://www.cfdonline.com/Forums/ope...ntensors.html http://www.cfdonline.com/Forums/ope...ncemodel.html Since no one in this forum could help me so far, I thought it is perhaps a good idea to ask someone, really familiar with OpenFOAM, like you. I am aware that you are short of time, but I would really appreciate if you could perhaps help me or if you know any other persons, who already have implemented new turbulence models in OpenFOAM, whom I could ask for help. Thank you very much! 
This is the first answer from hjasak
Quote:

Hello Mr Jasak,
thank you for your answer. I apologize, because I forgot the definition of b, here is another implementation, where it probably becomes clearer, what I did.: Code:
volSymmTensorField b = R_/(2.0*k_)1.0/3.0*I; (skew(fvc::grad(U_)))&b) = (b&skew(fvc::grad(U_))^T) I could not find any mathematical explanation for this. I think, there is no way to express Code:
+C5_*k_*(b&W+W&b) Concerning your summary 1: How can I write a source term such that I never use full tensors? I really appreciate your answering. Thank you very much! Many Greatings from California 
This is the second answer from hjasak:
Code:
Close + you should really know this. Here goes: 
Hi Hrv,
thank you for your answer. With your help, I could now finally solve my problem and compile my new turbulence model. Unfortunately, it doesnt work properly. When I start a calculation it looks good, but after some iterations the Courant number begins to "explode", that means it is geting bigger and bigger and then I always get an "floating point error". Me and my advisor think, it has something to do with how OpenFOAM calculates the velocity gradient at the wall (because the model is sensitive to that). I tried to search in the OpenFOAM source code, where the implementation of the near wall velocity gradient is but I could not find it. Do you know where this is implemented or how OpenFOAM calculates the near wall gradients? Thank you very much for your help! PS: I hope the weather in London is better than it is here (because it is raining in san francisco)! 
This is the third answer from hjassak:
Quote:

Well, I now discovered that the implementation of the last summand of my model (see equ. (10) in first attached file of first post) is still wrong. Now, here is the dilemma:
I am very greatful for every hint or comment, that can help me with this. Thanks a lot! 
I found the error. My translation from Einsteinian tensor notation into matrix notation was not entirely correct. Instead of
symm(b&W+W&b) it must be symm(b&W+W&b) because the kindex of the W(jk) in b(ik)W(jk) is written on the outer side, W must be transposed, so b(ik)W(jk) becomes b&W(transposed) instead of b&W. W is skewsymmetric and for these matrixes, the transposed is the same as the negatived (W(transposed)=W). Thus, b&W(transposed) becomes b&W. The whole thing (b&W+W&b) then is symmetric and by applying the symm() function on that, we get only the 6 upper triangle parts of this 9component matrix and are therefore able to use it within the fvSymmMatrix template. Thus, my problem is solved. Thanks a lot for all your help, and special thanks to hjasak! Have a nice weekend:) 
Hi again,
as stated above, my model is now working, which means that my calculations are actually converging now. Nevertheless the results do not yield what we expected. After running several checks, my Professor and me got the impression, that the bad results are probably caused by how OpenFOAM treats the near wall region. I tried to find out how OpenFOAM does this and I found something like Code:
// Calculate nearwall velocity gradient I tried to get an understanding of how this works (how it calculates the velocity gradient at the wall) but I somehow got stuck somewhere in the source Code behind that. So my question is: Does someone know how OpenFOAM calculates the velocity Gradient at the wall? I appreciate all your help, thank you very much! 
Hi Sven,
Do you use wall functions, or do you resolve also the viscous sublayer? The gradient at the wall is probably computed in the simplest possible way: value at wall, value in first cell centre, and distance between them, etc... This will, of course, give you a true wall gradient only if the velocity profile is linear, i.e. for y+ ~ 1. /Ola 
Hi Olwi,
I do use wall functions. The distance of my first cell to the wall is about y+~30 to bridge the viscous sublayer with the wall function. That means I do not resolve the viscous sublayer. If the velocity Gradient would be calculated by OpenFOAM as you say (by only using a simple difference), that would in my opinion be badly wrong. Can it be that the velocity gradient is somehow calculated by the wall function? 
Hi Sven,
You should definately look through the source code more, to understand the way the turbulence models already supplied deal with the wall region and wall functions. From the postings above I get the feeling you have spent a lot of time on the core flow, but not so much on the wall treatment. A wall function is an intrinsic part of a turbulence model. What do YOU want to do with the wall region? Does your model specify a certain wall function? How does the analytic velocity profile of the wall function enter you production term and pressurestrain terms? Do you integrate the expression for the velocity profile over the first cell for these terms, to account for the "average" happenings there? Sorry, I'm not helping much I'm afraid... ;) I've noticed that OF 1.5 and 1.6 deal with wall functions differently. The 1.6way is perhaps more "structured" and organized, but I think the more pragmatic way it's done in 1.5 is easier to understand, if you look through the code. Remember that wall friction tau=mu*dU/dn at the wall. In 1.5, at least, the wall function in the momentum equation is accounted for by modifying the value of mu (!) at the wall, rather than using the wall function to compute the gradient dU/dn differently! This is simply easier to implement, to avoid messing around inside the boundary condition classes (which is where the gradient is computed, I think). In 1.6 I noted that boundary conditions are involved, but my guess is that the principle is the same, but the trick has moved from one part of the library to another. Good luck with the sourcecode digging! /Ola 
Hi olwi,
first of all thanks for your answer. I implemented my own turbulence model by editing the existing LaunderGibsonRSTM from OpenFOAM 1.6. That means, I replaced the existing equation for the Reynolds Stresses with the new equation for my new turbulence model. That is basically the only change I made. For all the rest I used the source code given by the LaunderGibsonRSTM, that means I also inherited how this code treats the wall. I think the wall treatment is implemented like this (these are some lines form the LaunderGibsonRSTM.C file): Code:
// Correct wall shear stresses Thanks a lot! yours Sven 
All times are GMT 4. The time now is 00:22. 