[sven]

I want to implement a new turbulence model (Dafalias-Younis) into OpenFOAM. This turbulence model is a Reynolds Stress model, which models the pressure-strain 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!

[sven]

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:

 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_*(bb-1.0/3.0*trace_bb*I)
        +C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_)))&b)     
        
     );

Unfortunately I cannot compile this edited file anymore. I tried a lot of things and the problem seems to be, that the last term in the equation of the edited model:

 +C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_)))&b)

uses the skewsymmetric-part of the velocity gradient. Therefore this is a full matrix (9 components), while all the other components only consist of 6 components (since they are symmetric). I could not find a solution to solve this problem for weeks, although I already posted some threads in the forum.

http://www.cfd-online.com/Forums/ope...n-tensors.html

http://www.cfd-online.com/Forums/ope...nce-model.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!

[sven]

This is the first answer from hjasak

Relatively easy: since you are solving the equation for a symmetric tensor, all your source terms must be symmetric tensors. The second thing is that there is something badly wrong with your model, but let's go through the mechanics first:

+C5_*k_*(b&skew(fvc::grad(U_)))+(skew(fvc::grad(U_ )))&b)

Ypur problem is that
(b&skew(fvc::grad(U_))) and the other one are full tensors. In that case, recognise:

(skew(fvc::grad(U_)))&b) = (b&skew(fvc::grad(U_))^T)
and now you have
b & (skew(fvc::grad(U_)) + skew(fvc::grad(U_)))^T

The lot in brackets is now symmetric (2*symm(...)) and all is well.

However, your model is sick or you mis-understood it because symm(skew(grad U)) = 0.

Summary: find yourself a symmetric tensor in the rhs and write the source terms such that you never use full tensors.

Summary 2: re-visit the paper and find out what went wrong

Summary 3: grad(U) is a second-rank tensor, dotted with b gives a vector. My guess is that you will have a diadic (bb), but it's only a guess.

Good luck with your work; feel free to publish my response.

Hrv
__________________
Hrvoje Jasak

[sven]

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.:

  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_*(bb-1.0/3.0*trace_bb*I)
        +C5_*k_*(b&W+W&b)        
     );

Here, b is a second rank tensor and W as well. I don't understand your proposal:
(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

 +C5_*k_*(b&W+W&b)

by only using symmetric matrices. Is there any way to get OpenFOAM to deal with non-symmetric matrixes in this equation?
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

[sven]

This is the second answer from hjasak:

Close + you should really know this.  Here goes:

- a simple re-ordering:

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

[sven]

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)!

[sven]

This is the third answer from hjassak:

I would say you simply need to be more careful when running the code. It is known that RSTM is difficult to run and people (including myself) typically run an eddyviscosity model first, and then restart full RSTM with the eddy viscosity result as the initial guess.

You also need to be careful with under-relaxation factors etc: this will allow you to keep the cude running and get a good result.

Best,

Hrvoje
__________________
Hrvoje Jasak

[sven]

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:

1. First of all, this summand is written in Einsteinian notation (see again first attached jpg) and I need to transfer it in Matrix notation. I did translate the individual terms like this: b_ij is b (b=R/(2k)-1/3*I; S_ij is S (S=symm(fvc::grad(U_)) ); and W_ij is W (W=skew(fvc::grad(U_)) );
   The full term, I translated as:
   C5*k_*(b&W+W&b). I am not completely sure if this implementation is correct like this. Thus this is the first uncertainness.
2. If the translation from 1. is right, there occurs the next problem. Since, W is a non-symmetric tensor, so is b&W and W&b. The template in which the equation, and therefore also this term is included (see 2nd post), seems to be only working for <fvSymmTensorMatrix>. Since b&W and W&b is not symmetric, I am not able to compile the code with the newly included term. The problem is, that, except of b&W and W&b all the terms in the equation within the template are symmetric. OpenFOAM stores symmetric matrixes as 6 components (the other three components can be deduced from these) and unsymmetric matrices as 9 components. In my opinion that is why I can't compile the code.
3. I tried a lot of different things to convert symmetr. matrices in nonsymmetric ones and vice versa, but I did not succeed so far. I do not know if there is another way to force OpenFOAM to deal with non-symmetric tensors in this template.

So, these are mainly the problems I have while implementing this new turbulence model in OpenFOAM. At the moment I am kind of desperate, because I cannot solve this problem now for several weeks.
I am very greatful for every hint or comment, that can help me with this. Thanks a lot!

[sven]

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 k-index 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 skew-symmetric 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

[sven]

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

 // Calculate near-wall velocity gradient
                tensor gradUw
                    = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];

in the .C file of my turbulence model.

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!

[olwi]

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

[sven]

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?

[olwi]

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 pressure-strain 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.6-way 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 source-code digging!

/Ola

[sven]

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):

// 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 near-wall velocity gradient
                tensor gradUw
                    = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];

                // Calculate near-wall shear-stress 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();
            }
        }
    }

Unfortunately I do not completely understand this part of the source code, but it seems like they calculate a near wall velocity gradient, with which a new wall shear stress tensor is calculated. If you understand this part of the code better than I do, I would be happy if you could tell me what exactly is done here!
Thanks a lot!

yours Sven