
[Sponsors] 
September 17, 2009, 15:18 
Implementing new turbulence model

#1 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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! 

September 17, 2009, 15:20 

#2 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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 ( fvm::ddt(R_) + fvm::div(phi_, R_) // fvm::laplacian(Cs_*(k_/epsilon_)*R_, R_)  fvm::laplacian(DREff(), R_) + fvm::Sp(K1_*epsilon_/(2*k_)+K1Star_*G/(2*k_),R_) == P +(1.0/3.0)*(K1_*epsilon_+K1Star_*G)*I(2.0/3.0)*epsilon_*I +(C3_C3Star_*(pow(trace_bb,0.5)))*k_*S +C4_*k_*((bS+Sb)2.0/3.0*tr(bS)*I) +K2_*epsilon_*(bb1.0/3.0*trace_bb*I) +C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_)))&b) ); Code:
+C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_)))&b) Conversion of tensors Editing RSTM Turbulence Model 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! 

September 17, 2009, 15:22 

#3  
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
This is the first answer from hjasak
Quote:


September 17, 2009, 15:24 

#4 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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; volSymmTensorField bb = b&b; volScalarField trace_bb = tr(bb); volSymmTensorField S = symm(fvc::grad(U_)); volTensorField W = skew(fvc::grad(U_)); volSymmTensorField bS = b&S; volSymmTensorField Sb = S&b; tmp<fvSymmTensorMatrix> REqn ( fvm::ddt(R_) + fvm::div(phi_, R_) // fvm::laplacian(Cs_*(k_/epsilon_)*R_, R_)  fvm::laplacian(DREff(), R_) + fvm::Sp(K1_*epsilon_/(2*k_)+K1Star_*G/(2*k_),R_) == P +(1.0/3.0)*(K1_*epsilon_+K1Star_*G)*I(2.0/3.0)*epsilon_*I +(C3_C3Star_*(pow(trace_bb,0.5)))*k_*S +C4_*k_*((bS+Sb)2.0/3.0*tr(bS)*I) +K2_*epsilon_*(bb1.0/3.0*trace_bb*I) +C5_*k_*(b&W+W&b) ); (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 

September 17, 2009, 15:26 

#5 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
This is the second answer from hjasak:
Code:
Close + you should really know this. Here goes:  a simple reordering: b & W + W & b = b & (W + W^T); By definition:  (W + W^T) is a symmetric tensor.  b can be written as a sum of symmetric and antisymmetric part b = 1/2(b + b^T) + 1/2(b  b^T) Now the first part (symm(b)) is symmetric and the second (skew(b)) antisymmetric. For ALL tensors: a dot product between a symmetric and an antisymmetric tensor is zero. Thus: b & (W + W^T) = b & (2 symm(W)) b = symm(b) + skew(b) and skew(b) & (2 symm(W)) = 0!!! Thus: b&W + W&b = symm(b) & 2 symm(W) and you are on velvet! If you're going to San Francisco Be sure to wear some flowers in your hair If you're going to San Francisco You're gonna meet some gentle people there Enjoy, Hrv __________________ Hrvoje Jasak 

September 17, 2009, 15:26 

#6 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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)! 

September 17, 2009, 15:32 

#7  
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
This is the third answer from hjassak:
Quote:
Last edited by sven; September 17, 2009 at 16:50. 

September 17, 2009, 16:49 

#8 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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! 

September 18, 2009, 20:50 

#9 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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 

September 29, 2009, 21:09 

#10 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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 tensor gradUw = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei]; 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! 

October 1, 2009, 17:46 

#11 
Member
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8 
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 

October 2, 2009, 20:48 

#12 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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? 

October 5, 2009, 15:41 

#13 
Member
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8 
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 

October 12, 2009, 13:18 

#14 
Member
Sven Winkler
Join Date: May 2009
Posts: 70
Rep Power: 8 
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 forAll(patches, patchi) { const fvPatch& curPatch = patches[patchi]; if (typeid(curPatch) == typeid(wallFvPatch)) { symmTensorField& Rw = R_.boundaryField()[patchi]; const scalarField& nutw = nut_.boundaryField()[patchi]; vectorField snGradU = U_.boundaryField()[patchi].snGrad(); const vectorField& faceAreas = mesh_.Sf().boundaryField()[patchi]; const scalarField& magFaceAreas = mesh_.magSf().boundaryField()[patchi]; forAll(curPatch, facei) { // Calculate nearwall velocity gradient tensor gradUw = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei]; // Calculate nearwall shearstress tensor tensor tauw = nutw[facei]*2*symm(gradUw); // Reset the shear components of the stress tensor Rw[facei].xy() = tauw.xy(); Rw[facei].xz() = tauw.xz(); Rw[facei].yz() = tauw.yz(); } } } Thanks a lot! yours Sven 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
SimpleFoam case with SpalartAllmaras turbulence model implemented  nedved  OpenFOAM Running, Solving & CFD  2  November 30, 2014 23:43 
EulEul flow, kekpepTheta Turbulence model  us  FLUENT  5  April 5, 2011 02:29 
KOmega Turbulence model from wwwopenFOAMWikinet  philippose  OpenFOAM Running, Solving & CFD  30  August 4, 2010 10:26 
validating turbulence model on flow around a car  Pedro  CFX  1  February 20, 2008 17:32 
SSG Reynolds Turbulence Model  Georges  CFX  1  February 28, 2007 17:15 