[gu1]

Hello,

I am in need of help (OpenFOAM 5.0), unfortunately I do not know much where to start and so I am here in the forum. I need to implement a function in the k-Epsilon model, called the Peterlin function (attached image).
'L' is a constant and 'C' is a matrix. How could I write the code and tell OpenFOAM that I want the elements of the main diagonal (in the division)? And how will it interpret that the process should be in steps (first the C11 ... then the C22 ...)?

Understanding this is very important to me.

[usv001]

If C is a square matrix, then you can loop over the elements to get the diagonal as a field like this:

scalarField diag(C.n());

forAll(diag, k)
{
    diag[k] = C[k][k];
}

If C is a tensor, then you can create the diagonal field by calling the individual components like this:

scalarField diag(3);

diag[0] = C.xx(); diag[1] = C.yy(); diag[2] = C.zz();

USV

[gu1]

Originally Posted by usv001

If C is a tensor, then you can create the diagonal field by calling the individual components like this:

scalarField diag(3);

diag[0] = C.xx(); diag[1] = C.yy(); diag[2] = C.zz();

USV

Yes, it is a tensor. But I imagine I can call Ckk from tr(C). Do you agree?

I have one more question and I hope you can help me. Look this:

Link
// Stress transport equation
fvSymmTensorMatrix AEqn
(
fvm::ddt(A_)
+ fvm::div(phi(), A_)
==
twoSymm(C)
- fvm::Sp(varf_/lambda_, A_)
+ (a/lambda_) * Ist
);

This part of the code is similar to what I described in the image, but before equality I have a conformation tensor.

The point I want your attention is, this code will give me 6 values (since the matrix is symmetric), consequently fvMatrix == fvMatrix. Therefore, I can not enter a scalar value into this equation, right?!

I ask this because I have another question and maybe you can help me too. When I analyze the code:

//kEpsilon.C
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(k_)
+ fvm::div(k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_()/k_(), k_)
);

G is a scalar according to the link, right?

So... I'm inferring the logic I wrote above (fvMatrix == fvMatrix). Now I have a '' == scalar - fvMatrix_ ", or no?!

...and another detail is, when I call k_(), I'm calling a scalar or a matrix?
What does fvScalarMatrix do?

All this questioning is because I tried to make this modification in the code:

// Stress transport equation
fvSymmTensorMatrix AEqn
(
fvm::ddt(A_)
+ fvm::div(phi(), A_)
==
twoSymm(C)
+ gama2 * k_()/lambda_ //gama2 and lambda_ are constant.
- fvm::Sp(varf_/lambda_, A_)
+ (a/lambda_) * Ist
);

and OpenFOAM did not accept.
All these questions are very important to me, believe me... if you can help me, it will make a big difference.

[gtarang]

fvScalarMatrix constructs a matrix for the scalar equation 'k'. So in the 'k' equation, 'G' is scalar or Production Term (Sij.Sij, Sum of squares of the components of rate of strain tensor).
Coming to the second equation which you post, here you are constructing a solver matrix for Tensor variable, hence all the terms should be of Tensor type (dimensional balance)
Now coming to the modification you are trying to do for the k-Epsilon, I think you will be adding some terms to k equation from Cij, and you will get Cij from other equation. For this, I will suggest you to write the 'k-e' equation without any notation form and see the indexes of Cij coming in. The thing is 'k-e' are basically scalar transport equation, so the term which your adding should be scalar. Now, to extract the scalar terms from the tensor, have a look in the programming manual of Openfoam, you will find appropriate function.


-
gtarang

[usv001]

Let me try to answer your questions in a systematic manner to the best of my knowledge.

Yes, it is a tensor. But I imagine I can call Ckk from tr(C). Do you agree?

tr(C) gives you the sum of the diagonal components not the individual components.

The point I want your attention is, this code will give me 6 values (since the matrix is symmetric), consequently fvMatrix == fvMatrix. Therefore, I can not enter a scalar value into this equation, right?

Yes, this is fvSymmTensorMatrix. So, all the terms should be symmetric tensors too as pointed out by Tarang.

G is a scalar according to the link, right? So... I'm inferring the logic I wrote above (fvMatrix == fvMatrix). Now I have a '' == scalar - fvMatrix_ ", or no?!

Your thinking is right but, technically, G is a volScalarField. So, you have "== volScalarField - fvScalarMatrix". In the matrix problem, Ax = b, terms in G will be added to b.

when I call k_(), I'm calling a scalar or a matrix?

You're calling the currently stored value of k. So, it is a volScalarField too.

So, your modification,

// Stress transport equation
fvSymmTensorMatrix AEqn
(
      fvm::ddt(A_)
    + fvm::div(phi(), A_)
    ==
      twoSymm(C)
    + gama2 * k_()/lambda_ //gama2 and lambda_ are constant.
    - fvm::Sp(varf_/lambda_, A_)
    + (a/lambda_) * Ist
);

is incorrect since you have a volScalarField instead of volSymmTensorField. I am not sure about what equation your're trying to implement. Please post the equation you're trying to model.

USV

[gu1]

Hi usv001,

...thanks for your reply.

I want to insert the equation below. You can use this link as a guide.

The first terms of the equation are div(phi(), A_) and laplacian(nut(), A_). The rest I'm having trouble writing. As I said, the link serves as a reference.

I understand that to simplify the process, I could call Ckk from tr(C) and unfortunately I do not know if this would affect the results (maybe) but, the logic you showed me, I understand how complex to my knowledge of C++ and I do not know how to implement it. If you can help me.

So, Ckk is one of my biggest problems (...and this also involves the div and the laplacian) and also the second term after equality (whose are scalar values).. I was thinking of multiplying k_ * I, but, I dont know... or something like this: (scalar) * I.

EDIT1: How would you write the S² equation? volSymmTensorField twoD = twoSymm(fvc::grad(U()); ?

[gu1]

Hi usv001,

I've made some great progress but I still need some help because I do not know if what I'm doing is correct.

The solver compiles normally and when I ask to execute it it presents the floating point error, follows a piece of the log:

...
Pressure gradient source: uncorrected Ubar = 1.47222, pressure gradient = 0.292009
DILUPBiCGStab: Solving for Cxx, Initial residual = 1, Final residual = 4.79226e-10, No Iterations 3
DILUPBiCGStab: Solving for Cxy, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Cxz, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Cyy, Initial residual = 1, Final residual = 4.79226e-10, No Iterations 3
DILUPBiCGStab: Solving for Cyz, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Czz, Initial residual = 1, Final residual = 4.79226e-10, No Iterations 3
ExecutionTime = 18.94 s ClockTime = 19 s

Time = 19

DILUPBiCGStab: Solving for Ux, Initial residual = 0.0801566, Final residual = 8.68273e-07, No Iterations 2
DILUPBiCGStab: Solving for Uy, Initial residual = 0.0789024, Final residual = 8.57377e-07, No Iterations 2
DILUPBiCGStab: Solving for Uz, Initial residual = 0.0179611, Final residual = 3.17481e-06, No Iterations 1
Pressure gradient source: uncorrected Ubar = 1.47222, pressure gradient = 0.306078
GAMG: Solving for p, Initial residual = 0.260731, Final residual = 6.2926e-07, No Iterations 16
time step continuity errors : sum local = 2.81231e-10, global = -8.32353e-18, cumulative = 2.61429e-18
Pressure gradient source: uncorrected Ubar = 1.47222, pressure gradient = 0.298783
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5 Foam::fvMatrix<Foam::SymmTensor<double> >::solveSegregated(Foam::dictionary const&) at ??:?
#6 Foam::fvMatrix<Foam::SymmTensor<double> >::solve(Foam::dictionary const&) at ??:?
#7 Foam::fvMatrix<Foam::SymmTensor<double> >::solve() at ??:?
#8 Foam::FENE_P::correct() at ??:?
#9 ? at ??:?
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? at ??:?
floating point

I've checked my code and I believe that the floating point error may be happening due to the terms Cxy ... Cyz... anyway, these terms are zero due to an '' adaptation '' I made in the code ... but it's an addiction for lack of knowledge and this is where I need help. Here is the snippet of my code that refers the first 3 members to the equation shown in the previous comment.

dimensionedSymmTensor Ist
(
"Identity",
L2_.dimensions(),
symmTensor::I
);

volSymmTensorField S(symm(L));

volScalarField doubleS = S && S;

tmp<fvSymmTensorMatrix> CEqn
(
fvm::div(phi(), C_)
+ fvm::laplacian(nut, C_)
==
2. * (lambda_/sqr(fP_) * doubleS * Ist) //this is my "adaptation"
);
CEqn.ref().relax();
solve(CEqn);

The code only shows the values of Cxx, Cyy and Czz due to the term Ist (my adaptation). It was the only way I found to make a scalar (since the product of two symmetric matrices is a scalar) enter as a symmetric tensor in the equation.

I need to learn how to make the product of the tensor strain rate (Sij), because in the equation it is squared.
I believe that if I do this correctly, the values Cxy ... Cyz ... will reappear in the log.

[usv001]

Hi Guilherme,

Just to clarify, the link you provided has a 'fvSymmTensorMatrix AEqn' in which all 6 components are solved. In the equation you provided, the variable to be solved is mentioned as which should by right refer to the trace of C which is a scalar. However, you seem to want to solve for the diagonal components individually. Could you please clarify this issue? If you only need to solve for 3 components (or even 1), you can create separate fvScalarMatrix for each component and solve them individually before combining the solution into the volSymmTensorField. That would surely avoid the off-diagonal components.

USV

[gu1]

Originally Posted by usv001

However, you seem to want to solve for the diagonal components individually. Could you please clarify this issue?

Yes, that is my goal. Unfortunately I do not know how to solve it.
If you could help me, I would be very happy.

[usv001]

You could try something like the following to solve for individual diagonal components:

// labelList corresponding to diagonal components in symmTensor
labelList cmpts(3);
cmpts[0] = 0; // xx
cmpts[1] = 3; // yy
cmpts[2] = 5; // zz

forAll(cmpts, i)
{
    // extract individual component
    volScalarField Ci
    (
        C_.name() + name(cmpts[i]), // rename 
        C_.component(cmpts[i])
    );

    // solve for individual component
    fvScalarMatrix
    (
        // your eqn for Ci
    );

    // replace back with solved component
    C_.replace(cmpts[i], Ci);
}

However, I still don't understand what happens to the off-diagonal components since you're not solving for them...

USV

[gu1]

However, I still don't understand what happens to the off-diagonal components since you're not solving for them...
USV

Hello UNV, thank you for your response.
Apparently the off-diagonal are not solved, according to the equations attached below.

I understand that by using Ckk it is only referring to the elements of the main diagonal.

Taking advantage of it, is it possible to insert these equations directly into the solver simpleFoam and thus, to simultaneously use the turbulence models already implemented in the OF?

[gu1]

If I want to calculate a dumping function that is a function of uTau, how do I insert the uTau equation into the kEpsilon.C file??

[usv001]

Hello Guilherme,

Originally Posted by gu1

Taking advantage of it, is it possible to insert these equations directly into the solver simpleFoam and thus, to simultaneously use the turbulence models already implemented in the OF?

I see no problem doing that.

USV

[gu1]

Originally Posted by usv001

Hello Guilherme,

I see no problem doing that.

USV

USV,

I wrote this snippet of code:

    Dxy_
    (
        IOobject
        (
            IOobject::groupName("Dxy", U.group()),
            runTime_.timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        //symm(fvc::grad(U))().component(symmTensor::XY)
        U.mesh(),
	dimensionSet(0, 0, -1, 0, 0, 0, 0)
    ),

    volSymmTensorField D = symm(gU);    
    
    forAll(D, i)
    {
    	Dxy_[i] = D[i].xy();
    }

...and I succeeded in what I wanted.
Unfortunately, however, the gradient of xy on the wall is giving zero and this causes the solver to give a floating-point error.

log: LINK

Could someone help me with how to solve this problem?

Best.

*By the way, if I want to enter the value of Dxy in a tau (symmetric) equation, how would you recommend that I do?

[usv001]

Hello Guilherme,

Your code only seems to be transferring the internalField only. To transfer the boundary values, try the following:

forAll(D.boundaryField(), patchi)
{
    Dxy_.boundaryField()[patchi]
        == D.boundaryField()[patchi].component(symmTensor::XY);
}

This should transfer the boundary values and, hopefully, that solves your problem.

Please post the equation for tau in code form so that people know what exactly you are trying to do because Dxy_ (a volScalarField) could be introduced many ways, as a source term, as a coefficient multiplying another term, etc. in an fvSymmTensorMatrix.

USV

[gu1]

Originally Posted by usv001

Please post the equation for tau in code form so that people know what exactly you are trying to do because Dxy_ (a volScalarField) could be introduced many ways, as a source term, as a coefficient multiplying another term, etc. in an fvSymmTensorMatrix.

USV

Hi,

I've tried writing the code that sent me up this way:

    volSymmTensorField D = symm(gU);    
    
    forAll(D, i)
    {
    	Dxy_[i] = D[i].xy();

    	forAll(D.boundaryField(), patchi)
        {
            Dxy_.boundaryField()[patchi]
                == D.boundaryField()[patchi].component(symmTensor::XY);
        }
    }

but... :

error: no match for ‘operator==’ (operand types are ‘const Foam::fvPatchField<double>’ and ‘Foam::tmp<Foam::Field<double> >’)
                 == D.boundaryField()[patchi].component(symmTensor::XY);
                 ^

Speaking about my code I have attached an article that I am referring to below, I believe it will be easier for you to visualize the equations.

LINK

I'm trying to write equations 8 and 11 described in this article. These are the set of the conformation tensor (equation 18) that compose the tau equation (equation 19). The way I've written both the last two is below.

    volScalarField fP = (L_-3.)/(L_-tr(C_));

    // Stress transport equation
    tmp<fvSymmTensorMatrix> CEqn
    (
        fvm::div(phi_, C_)
     ==
        M
      + NLT
      - fvm::Sp(fP/lambda_, C_)
      + (3./lambda_) * I
    );
    CEqn.ref().relax();
    solve(CEqn);  
    
    tau_ = (nuP_/lambda_) * (fP*C_ - I);    
    tau_.correctBoundaryConditions();

My difficulty is in modeling equation 8 (M), for lack of knowledge and necessary guidance.

It seems to me that M is related to: volSymmTensorField D = symm(fvc::grad(U()));
...but I dont know.

Another way in which it is possible to write equation 8 is one that was presented by a precursor author to the study presented. I attached the equation as an image.
So, for not knowing how to write in a way (eq. 8), I'm limiting myself to writing the problem in a 2D way, using the attached equation, it uses the squared gradUxy (this is possible since it's a scalar). And that's why I asked you for help.
This leads me to question that if I make this adaptation, the equation of the symmTensor does not necessarily need to be fvSymmTensorMatrix, it would only be fvScalarMatrix, since it will solve only the xy component and consequently it will only give me tauxy (which is what I need for a problem 2D).
I want to remember that NLT is in function of M and both are not demonstrated in the code that I attached because not yet how to write them.
*Incidentally, equations 22 and 24 described in the article are also dependent on M, but in this case, I BELIEVE that they are only dependent on the xy component of the U gradient.

For me, the ideal would be to write the equation 8 and thus, to be able to model a 3D problem.

All this argument is in the hope that you can help me.

PS: Incidentally, if I can ask for more detail, in case I want to write gradU², that is, sqr(fvc::grad(U())), how would it be done?