Access and Modify Turbulence Viscosity 'nut'

[JihooKang-KOR]

Dear Andre,

I would like to report the current status of the 'nut' issue. I thought of some possibilities to access and modify the whole field of turbulence viscosity as follows:

1. Currently, in the new wall function code 'nutUIntegralSpaldingWallFunctionFvPatchScalarField.C', we can access the whole 'nut' field by tmp<volScalarField>& tnut = turbModel.nut(), but the modification of the value except the wall region is impossible because it is a read-only mode. Therefore, I tried to find where the variable 'nut' is located (in No.2).

2. The class 'nutUIntegralSpaldingWallFunctionFvPatchScalarField' is inherited from the class 'nutWallFunctionFvPatchScalarField'. In the class 'nutWallFunctionFvPatchScalarField', there is the function called 'nutw', and it referred to 'turbulenceModel::nut()' by using turbModel.nut()().boundaryField()[patchi]. In the following link, we can see the diagram of caller graph for 'nut'. https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1turbulenceModel.html#a172707e93d70ead6d609ebf6138e46d9 Among them, I chose one method called 'Foam::RASModels::continuousGasKEpsilon::correctNut' because the name of the method seems to be related to accessing the 'nut' field. In this method (the file 'continuousGasKEpsilon.C'), there is nothing special, but only the access of 'liquidTurbulence.nut()'. Thus, I checked the same method ('correctNut()') in the parent class 'kEpsilon'. I could find that there is an access of 'this->nut'. However, there is no variable of 'this->nut' in the class, and therefore I checked the parent class 'eddyViscosity'. Finally, 'volScalarField nut' was found in 'eddyViscosity.H' as a protected variable in the class. In conclusion, we need to inherit the class 'eddyViscosity' in order to access the 'nut' variable. From this point, I couldn't decide how to implement it since if we inherit 'eddyViscosity' as well as 'nutWallFunctionFvPatchScalarField' simultaneously in the class 'nutUIntegralSpaldingWallFunctionFvPatchScalarField', there might be problems for encapsulation and compilation of the code.

3. In 'nutWallFunctionFvPatchScalarField.C', there is the method 'updateCoeffs()'. In the method, the code operator==(calcNut()); can be found. I thought we would be able to change the type of 'operator==' from passing 'scalarField' to passing 'volScalarField', but I lost the track to find the base operator overloading method. In addition, even if I find where it is, it would not be easy to update the whole 'nut' field only by changing 'operator=='.

I couldn't find any clear ways to modify 'nut' directly due to some limitations, but we might give it a try to write a code in terms of No.2 or No.3 if the approach is reasonable enough.

Thank you for reading this comment.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for digging into the BC code. I checked the following three options to implement the wall model:

• standard wall function as you discussed in your post
• function object
• finite volume option

I came to the following conclusions, which agree with your observations:

• the wall function structure/fixedValue BC is not suitable for two reasons: 1) it is designed to correct only patch values 2) we cannot easily modify the convective flux and we cannot ensure that the correction is applied at the right time to the correct field
• function objects are designed to be for post processing
• fvOptions is probably the right infrastructure for the implementation, since we can manipulate the matrix coefficients at the right time in the right place

Implementing a new turbulence model would be the right way to go if we only had to modify nut. However, we also need to modify convective fluxes. I do not have any experience yet with implementing new fvOptions, so I suggest hardcoding the corrections into the solver for now:

• create a new solver based on simpleFoam; something like wmSimpleFoam (wall modeled simpleFoam)
• include the source code you have already written for the boundary condition to compute uTau
• refer to my PhD repo on how to implement the convective and diffusive corrections for species transport; you'll need to do the same for the UEqn in simpleFoam; here is the tricky part that I didn't yet have to deal with: the nut computation is hidden inside the turbulence model, so you have to analyze turbulence->divDevReff(U) to see how to modify nut

That's all I have so far. Please, have a look at fvOptions, too. Best, Andre

[AndreWeiner]

Hi Jihoo, writing the former reply, in particular, the last bullet point, has given me second thoughts. Probably, we should implement a new incompressible::turbulenceModel, where the correction of nut is done in turbulence->divDevReff(U) and another function turbulence->divEff(U) replaces the fvm::div(phi,U) term. So, the idea would be to implement a new type of incompressible turbulence model. Best, Andre

[JihooKang-KOR]

Dear Andre,

Currently, I am working on implementing the nut calculation part to the new turbulence model. I committed once for this, but there are some problems as follows:

1. The code is moved from the wall function nutUIntegralSpaldingWallFunction to the new model wmSpalartAllmaras, and I tried to update the value at the exact wall patch. However, this->nut_ has only the internal field which can be manipulated. (There is no boundary field in the nut_ variable itself.) We can access the wall patch by this->nut_.boundaryField()[patchi], but we cannot manipulate these values because they are accessed as constant. Therefore, calcNut() in the new model calculates nut at the first cell center normal to the wall, and pass it to this->nut_.correctBoundaryConditions() to update the wall patch. Due to this problem, the Cf values are smaller than the new wall function as follows (orange line): The values were updated to the wall patch only when the boundary condition zeroGradient in 0/nut was applied. (If calculated is used, the values at the wall behave like no wall function.) It is not exactly the same with using wall function, but the BC is affected by the neighbor cell center values. For comparison, this is the Cf graph when the wall function nutUIntegralSpaldingWallFunction is used: Compared to this, the values of the first cell center based nut calculation are smaller.

2. I assumed that the boundary condition is updated by the code line this->nut_.correctBoundaryConditions() after the calculation in calcNut() since without this code, the values at the wall behave like no wall function. However, I couldn't find the original location of correctBoundaryConditions() anywhere. Now I don't know how to update the value at the wall using the value at the first cell center (interpolation or others?), and thus I need to find the location of the function.

3. I have a question about folder names. I set the folder name for solver and turbulence model as solver and turbulence_model because we currently have only one solver and one turbulence model. However, I'm not sure if we use more solvers or turbulence models. In this case, I would like to know if I should change the folder names to solvers and turbulence_models.

These are what I found so far that I would like to report. After solving the problems, I can proceed with finding the correction method for fluxes.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear Andre,

It seems that the problem No.1 mentioned above might be solved with the commit 304aaa7. By employing a reference variable (with &) and the function boundaryFieldRef(), I've got the similar results to the new wall function as follows.

Therefore, No.2 (about the location of correctBoundaryConditions() function) and No.3 (about the folder name in the repository) are remaining to be dealt with.

From now on, I will start to learn more about how SIMPLE algorithm works and the method of flux correction.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for the update.

2) the idea is to have the turbulence model compute and set nut at the wall patch; then the calculate BC (doing nothing) should be the right choice, and you don't have to worry about when or where correctBoundaryConditions() is called. 3) the folder names a fine as they are for now

If I understand it correctly, as of now you only correct nut at the wall patch. The next step would be to modify the nut value also at the first cell face normal to the wall patch. If that works, we move on to correcting the convective fluxes (we also need a new test case).

Best, Andre

[AndreWeiner]

Hi @JihooKang-KOR, here are some thoughts on how to resolve the implementation issues. I suggest a top-to-bottom approach (starting at the top level of the solver and then working towards implementation details). So here is a rough outline:

• we implement the correction inside the solver (e.g. wmSimpleFoam)
• instead of calling turbulence->divDevReff(U) in UEqn.H we include the function's implementation at the same place inside the solver

  
  // from linearViscousStress.H

• fvc::div((this->alpha*this->rhothis->nuEff())dev2(T(fvc::grad(U))))
  • fvm::laplacian(this->alpha*this->rho*this->nuEff(), U)

• how to replace the term:
  • in createFields.H, create a surfaceScalarField called nuEff
  • before creating and solving the momentum equation, set nuEff to fvc::interpolate(turbulence->nuEff())
  • overwrite the face values of nuEff at the boundary patch and the first faces normal to the wall using the integral fitting/correction approach
  • the new form of the expression above should be

• fvc::div(nuEff, dev2(T(fvc::grad(U)))) - fvm::laplacian(nuEff, U)

  The resulting momentum equation should be:

  tmp tUEqn ( fvm::div(phi, U)

  • MRF.DDt(U)

  • fvc::div(nuEff, dev2(T(fvc::grad(U)))) - fvm::laplacian(nuEff, U)

    fvOptions(U) );

Here are some useful links:

• RAS turbulence model guide
• OF Programmer's guide; see p. 35 for the definition of various operators (div, laplace, etc.)
• linearViscousStress
• turbulenceModel

I suggest implementing these changes step by step, e.g., first, replace the term in the momentum equation without doing any additional corrections at and close to the wall; check that the results are unchanged compared to the default implementation; then add the corrections. There might be also implementation details that are still missing. We'll discuss them as they come up.

This is the trickiest part of the project, but I am confident that we'll get it done :-) Let me if you get stuck!

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

While modifying the momentum equation in UEqn.H, I faced a problem as follows.

First of all, the problem is caused by the term fvc::div. When I tried fvc::div(nuEff, dev2(T(fvc::grad(U)))) mentioned above, the error was given below:

error: no match for ‘operator-’ (operand types are ‘Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >’ and 
‘Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >’)

I think this was because of wrong number of parameters, and in turbulence->divDevReff(U) function, the equation is expressed differently. Therefore, I tried fvc::div(nuEff * dev2(T(fvc::grad(U)))) this time, and the error was given below:

error: no match for ‘operator*’ (operand types are ‘Foam::surfaceScalarField’ {aka ‘Foam::GeometricField<double, 
Foam::fvsPatchField, Foam::surfaceMesh>’} and ‘Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, 
Foam::fvPatchField, Foam::volMesh> >’)

It seems that we cannot use the * (multiplication) operator for surfaceScalarField in fvc::div(). Hence, I consequently tried the momentum equation as follows:

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      - fvc::div(turbulence->nuEff() * dev2(T(fvc::grad(U)))) 
      - fvm::laplacian(turbulence->nuEff(), U)
     ==
        fvOptions(U)
    );

Then, it was complied well, and the simulation result was the same for y+ = 100. But turbulence->nuEff() is the function for cell center values, so it seems that we might not be able to directly use surfaceScalarField in the momentum equation.

That is why I committed 0fa160c and 3cde2b3 to apply the corrected face values to the cell centers inside the turbulence model that we've created. If we cannot directly apply surfaceScalarField in the momentum equation in the new solver, we need to adapt the situation and use cell center values as OF originally uses. The result plot is as follows: There are some discrepancies of Cf for each y+, but no more giggling exists and the graph patterns are similar. In the notebook (the commit 3cde2b3), a brief explanation is available.

I'm not sure if we have to modify the file that contains the code of fvm and fvc to use surface values, but I would like to report this problem as soon as possible in the Github issue.

Thank you very much for reading this comment.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear @AndreWeiner,

I did two simulations as follows, but the results are unfortunately not what we expected.

1. The following momentum equation where surfaceScalarField nuEff = fvc::interpolate(turbulence->nuEff()).

   tmp<fvVectorMatrix> tUEqn
   (
   fvm::div(phi, U)
   + MRF.DDt(U)
   - fvc::div(turbulence->nuEff() * dev2(T(fvc::grad(U)))) 
   - fvm::laplacian(nuEff, U)
   ==
   fvOptions(U)
   );

   In this case, the original turbulence->nuEff() that is volScalarField is used in fvc::div, and the interpolated surfaceScalarField nuEff is used only in the laplacian term.

2. The momentum equation without the fvc::div term.

   tmp<fvVectorMatrix> tUEqn
   (
   fvm::div(phi, U)
   + MRF.DDt(U)  
   - fvm::laplacian(nuEff, U)
   ==
   fvOptions(U)
   );

Both simulations were done with fixedValue 0.0 for the bottomWall patch in nut.

The results of Cf are as follows (for y+ = 0.05 and 30 here):

As seen in the figure, the simulation No.1 fluctuates a lot, and the simulation No.2 cannot capture the original behavior (SA model without wall functions) of Cf.

Regarding residuals, I would like to show the example residual for y+ = 0.05 as follows: The case with the fvc::div term. The residual did not converge, but it soared and yielded error.

The case without the term. The residual converged, but we found in the Cf graph that the behavior couldn't capture the real one.

According to the results above, we need to come up with other ideas.

Meanwhile, I would like to ask one thing about the principle of the nut correction. I understood how it works for one face (e.g. wall patch), but I don't understand why the behavior of Cf will almost be identical if we correct nut not only for the wall, but also for the first cell face opposite to the wall. I thought Cf would differ from each y+ because u_tau of integral fitting method is different from the simulated u_tau when the height y exceeds a certain value (making u_int = u_sim by Secant method). Therefore, I think that correcting first cell face nut makes the result more robust since a mathematical calculation is involved in the nut field, but Cf might be different among y+ due to the different u_tau. Thus, if you could give me an insight for this question, I would appreciate it.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for the update. I would like to understand where the difference arises if we pass a surfaceScalarField to fvm::laplacian instead of a volScalarField. I expected the operator to do a simple fvc::interpolate of the nut field. Could you explore the actual computation a little bit more in detail? If you pass nuEff as volScalarField, do you get the same results as without replacing the terms (the original implementation)?

Regarding the last point, the idea is that we can reconstruct the correct uTau and Cf given only the correct cell-centered velocity (representing the average velocity in the cell) and the cell-width by making some assumptions about the velocity profile (Spalding's function). When starting the simulation, we do not know the correct cell-centered velocity, but we have to find it iteratively by solving the pressure and momentum equation. In the discretization of the momentum equation, we will introduce a significant interpolation error if we use simple linear interpolation. The main error sources arise from the terms containing the surface normal derivative of the velocity component tangential to the wall. Therefore, we aim to improve these two interpolations (at the wall and at the first cell face normal to the wall) using Spalding's function. In principle, you could also try to correct the second/third/fourth and so on cell face normal to the wall, but usually, correcting the first cell face normal to the wall is enough. It is also less demanding considering the search operations and the requirements on the mesh (having multiple prismatic cell layers). Without these corrections, the solution of the momentum equation would not yield a good approximation of the cell-centered velocity and we would not be able to reconstruct uTau/Cf. E.g., you could consider using Spalding's function only as a post-processing step. However, this would give you mesh-dependent results, which are not as good as they could be because the cell-centered velocity you pass to the correction has a large, mesh-dependent error.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I tested two settings as follows:

1. withfvc : Apply turbulence->nuEff() (which is volScalarField) to both of the div and laplacian terms.

   tmp<fvVectorMatrix> tUEqn
   (
   fvm::div(phi, U)
   + MRF.DDt(U)
   - fvc::div(turbulence->nuEff() * dev2(T(fvc::grad(U)))) 
   - fvm::laplacian(turbulence->nuEff(), U)
   ==
   fvOptions(U)
   );

2. nofvc : Apply turbulence->nuEff() only to the laplacian term. (The fvc::div term is commented out.)

   tmp<fvVectorMatrix> tUEqn
   (
   fvm::div(phi, U)
   + MRF.DDt(U)
   //- fvc::div(turbulence->nuEff() * dev2(T(fvc::grad(U)))) (Commented out)
   - fvm::laplacian(turbulence->nuEff(), U)
   ==
   fvOptions(U)
   );

Compared to the original Spalart-Allmaras model, they tend to behave the same except for y+ = 5 and 10. Even for y+ = 5 and 10, the behavior is not dramatically different. (The identical behavior when y+ = 0.05)

(The different behavior when y+ = 5 and 10)

The residuals behave almost the same, but the no fvc term case seems to be more stable and converges a bit faster. I couldn't find any distinguishable differences among them.

Consequently, we can conclude that directly inserting the div and laplacian terms to UEqn.H let Cf values behave almost the same (not 100% same) when we pass nuEff as volScalarField, whereas the cases of nuEff which is the interpolated surfaceScalarField couldn't properly capture the behavior of the original SA model. According to the result, I think OpenFOAM calculates the results properly when we pass volScalarField to the momentum equation. However, as you explained before, we will lose some information of nut when the corrected nut at faces is interpolated to cell centers. In this situation, I'm not sure if we can use more sophisticated interpolation schemes to keep the nut information of faces or we still keep thinking how to directly pass surfaceScalarField for nuEff to OpenFOAM.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, I think there is something wrong with your implementation or the compilation. I compared the original simpleFoam results with the modified momentum equation below and got exactly the same results for yp=30:

    // Momentum predictor

    MRF.correctBoundaryVelocity(U);
    surfaceScalarField nuEff = fvc::interpolate(turbulence->nuEff());
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      - fvm::laplacian(nuEff, U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
...

Did you maybe place the nuEff computation outside the SIMPLE loop? Can you please upload the last modification including the test case and the evaluation?

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

Thank you very much for your comment. I realized that I put the surfaceScalarField nuEff outside of the SIMPLE loop. I will put this inside of the loop as you implemented above, and report the result by tomorrow.

Regarding the last modification, I uploaded two test cases (with fvc term and without fvc term using volScalarField turbulence->nuEff()) in the cloud storage. (I guess you would like to compare two settings of volScalarField and surfaceScalarField as an input to the laplacian term. Therefore, I uploaded these cases to the cloud.) The new solver codes are inside each folder. In addition to that, I made one commit to upload notebooks about Cf and residuals for those two cases.

Thank you again for the correction.

Best regards, Jihoo Kang

[AndreWeiner]

Great, thanks! I'll have a look at the folders. Best, Andre

[AndreWeiner]

Hi @JihooKang-KOR, here is a short summary of the next steps we discussed today and how to implement them. On the top level, at each time step/iteration, we:

1. get the effective viscosity from the turbulence model and interpolate it to the cell faces (surfaceVectorField nuEff = fvc::interpolate(turbulence->nuEff());)
2. overwrite the entries in nuEff corresponding to the boundary faces and the first faces normal to the boundary faces
3. pass nuEff as fvm::laplacian(nuEff,U) (forget about the fvc::div(...) term for now)

The actual correction is computed and applied in the second step as follows:

0. the portion of the code to determine cell and face indices of adjacent cells should be placed before the time loop if possible, since they are do not change at runtime
1. in each time step/iteration 1.1 fit Spalding's law to the velocity value in each cell at the boundary (find uTau for a given cell-centered velocity magnitude and the cell width) 1.2 determine the slope of the velocity profile at the boundary and at the first face normal to the boundary; if x and y are the coordinates along and normal to the wall, respectively, Spalding's law is a relation y=f(U_x) for a given uTau; what we need for the correction is dU_x/dy|_y=0 and dU_x/dy|_y=yf (yf denotes the distance between wall face and the first face in the normal direction); we can approximate these derivatives by using finite differences; first, define a small distance much smaller than yf, e.g. scalar eps = yf/1000.0;; the formulas to compute the derivatives are dU_x/dy|_y=0 = (U_x(y=0)-U_x(y=eps))/esp and dU_x/dy|_y=yf = (U_x(y=yf-eps)-U_x(y=yf+eps))/(2esp); note that this definition yields negative slopes; to obtain U at a given distance, the relation y=f(U_x) has to be inverted numerically, e.g., Newton/Bisecion method 1.3 compute and overwrite the effective viscosity for all boundary faces and the corresponding faces in normal direction; the formula is nuEff = nu {mag(dU_x/d_y)}_spalding / (small + {mag(dU_x/dy)}_linear); nu denotes the molecular viscosity (turbulence->nu()); instead of computing dU_x/d_y with the standard interpolation, it should be also fine to use mag(fvc::snGrad(U)) since dU_x/d_y dominates the other contributions
2. create functions that write all relevant information to simple CSV files (example); if you switch to parallel execution, there are two working approaches - you can have each process write its own file, e.g., to processor0/100/wall_model.csv, and gather the data during post-processing, or you can gather the field on the master and write the gathered field to postProcessing/wallModel/100/wall_model.csv (partial example); either way is fine - implement whatever is easiest for you; relevant information might be the new effective viscosities, the original values of nuEff, the cell-centered Velocity, the distance of the first cell face and the cell center, the velocity surface normal gradient obtained by linear interpolation, the coordinates of the boundary faces and the first face in normal direction, uTau, ...; you may also find some useful snippets in this file

Let me know if I forgot something or if questions/problems come up!

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I revised the nuEff correction part as per the above steps and the code was compiled well. However, after the temporary simulation with the condition of y+ = 100 and 500 iteration steps, I found that the simulation didn't go well as follows:

1. I calculated tauw with the formula nuEff_wall * magGradU or uTau^2, but both cases yielded the following figures. (In the code, I set nuEff as max(nu, calculation), and therefore if the value is less than the molecular viscosity nu, it becomes nu. In the log file, we can see nuEffface200 = 1.388e-05 for the case of using max function.) I think there wouldn't be proper results even if 5000 iterations were executed because the fluctuation appears a lot. <max function is applied> <max function is not applied>

2. I printed nuEff of the first face at the middle of the plate (nuEffface200) in a log file, but I found that this value went to the value that is smaller than the molecular viscosity nu after some iterations. I guess that it means nut becomes negative after some iterations, but this is different from the result of the new turbulence model which was implemented before.

I'm trying to find the reason, but I need more time to determine what is the problem. I upload the log files of the simulations and the result files for both versions ('surfaceFields.csv', by using very primitive printing code) here for reference.

log.wmSimpleFoam_yp100_1e-5_500iter.txt

surfaceFields.csv

log.wmSimpleFoam_yp100_1e-5_500iter_nomax.txt

surfaceFields_nomax.csv

<08.08.2021> In addition, I would like to add the result of 5000 iterations. <max function is applied> <max function is not applied>

As seen in the figures, they diverge before getting 5000 iterations. In the following residual figures, we can see that p field does not converge at all. <max function is applied for residuals> <max function is not applied for residuals> This might mean that we have to see the pressure field as well.

When I saw the log files, another problem appeared as follows:

3. For the case of using max function in the correction part, the velocity in x-direction was negative at the first cell center right before the core dumping.

4. For both cases, magGradU at the wall is smaller than that at the first face unlike the case of using the new turbulence model. I think the value at the wall should be larger.

I upload the log files again here for reference.

max_log.wmSimpleFoam_y100_5000iter.txt

max_surfaceFields_5000iter.csv

nomax_log.wmSimpleFoam_yp100_5000iter.txt

nomax_surfaceFields_5000iter.csv

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo,

I couldn't spot any obvious mistakes in the implementation. Here is what I suggest:

• could you please upload the test case you are currently using to do an individual serial test run?
• write out the surface fields at every n iterations; it is important to assess how the fields evolve in time
• write out the results of iterative processes and exception handling, e.g., if you limit gradients, write out the limited and unlimited values
• the effective viscosity should not be limited by the molecular one because the stresses might be sometimes over and sometimes underestimated; therefore, the effective viscosity might be larger or smaller than the molecular one; note that the effective viscosity is by construction larger or equal to zero since nu is larger than zero.
• start by correcting only the viscosity at the wall; only when that works properly, we also correct at the first cell face; you already made the wall-only correction work with the wall-function implementation, so you know what to expect - friction coefficients with the same trend as the normal wall function but slightly higher values
• make sure that the iterative algorithm for finding uTau always has plausible inputs and that it converges; especially the first iteration might pose a problem since the velocity field is initialized to the same value everywhere; one simple solution in that case would be to start with the correction only after the first x normal iterations; also, 10 integral points might be too little for yPlus=100

That's all that comes to my mind for now. We are close! Let's iron out the last couple of problems! Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I need more time to think about how to write out the related values because I couldn't find any rules of evolving the pressure and velocity. They fluctuate just randomly as I uploaded the log files before. However, I would like to report what I found so far in this issue as follows:

0. I changed the number of integral points from 10 to 50. uTau calculation yielded a proper value. The effective viscosity is not limited anymore.

1. First of all, only changing the wall patch yields the proper result as the new turbulence model, but when the first cell face correction is involved, the result starts to be fluctuated. This is caused by the fact that both of U fields U_face[patchFaceIDs[faceI]] and U.boundaryField()[surfaceID][faceI] refer to the same area (the wall patch), while they refer to the different areas (the first face or the first cell center) for the first cell face correction.

2. But I found that the calculation was properly executed (residuals converged well) when I made the code as follows:

   surfaceScalarField nuEff = fvc::interpolate(turbulence->nuEff());
   #include "integralFittingMethod.H"
   volScalarField nuEff_avg = fvc::average(nuEff);
   tmp<fvVectorMatrix> tUEqn
   (
   fvm::div(phi, U)
   + MRF.DDt(U)
   - fvm::laplacian(nuEff_avg, U)
   ==
   fvOptions(U)
   );

   I added a volScalarField nuEff_avg from nuEff by using fvc::average. (This is only the method that I know to change from surfaceScalarField to volScalarField.) Therefore, I could conclude that the momentum equation yields proper results when we pass all the cell center based fields. Here is the graph of the simulation that the above code was used:

As seen in the figure, Cf behaves a bit better than Cf of the new turbulence model except y+ = 5 and 10, but there might be a discrepancy between this result and what you expect since the average function was used. At least, I could find the implementation of the solver code is ok, but there would be problems in the way of the calculation in the momentum equation itself.

3. In conjunction with No.2, I checked how many number of elements are contained in fvm::laplacian(nuEff, U) and fvm::laplacian(nuEff_avg, U). Both of cases have 209825 elements. For the latter equation, it is obvious that two arguments are all volume fields. However, for the former equation, nuEff has 418720 elements because it is a surfaceScalarField. Thus, I guess that somewhere in OpenFOAM calculates or interpolates this surface field and apply to the laplacian function. But the simulation doesn't work in this case, whereas the simulation using volScalarField nuEff_avg = fvc::average(nuEff) converges well. For the case of inputting nuEff directly as fvm::laplacian(nuEff, U), we might need to think of dealing with x-direction faces as well (not the normal direction of the wall) to check the process of interpolation of face values.

4. For the case of using nuEff_avg, there was a discrepancy of the result, and therefore I tried to correct the second cell face as well. However, the result was not that good (too low Cf) compared to correcting until the first face. This was not intuitive for me, so I would like to inform of it.

I didn't write out the data except the above well worked case because I wanted to check the problem swiftly (did quick repetition of ./Allclean and ./Allrun). If you want to see any visible results of some cases, please let me know. Then, I will simulate and write out the data to csv files.

Thank you for reading this comment.

Best regards, Jihoo Kang

[AndreWeiner]

Dear Jihoo,

thanks for the update! Some comments on your observations:

1. the only velocity value you really need in the flat plate example is the velocity at the first cell center normal to the wall
2. if you apply fvc::average before fvm::laplacian, you average the nuEff surface field to the cell centers and then interpolate it back to the faces in the Laplace operator; the combination of averaging and linear interpolation simply smoothes out the nuEff field
3. by solving the momentum equation, you update the cell-centered values of the velocity field; in the finite volume discretization, the Laplace operator laplace::(nu,U) in a control volume is approximated as the sum over all faces composing the cell of nu on the face times the difference of the neighboring cell-centered values times some distance/area term (see eq. 1.34 of the OF technology primer); therefore, the face values of nu and not the cell-centered values enter into the coefficient matrix of the velocity equation; since the coefficients are computed as a sum of contributions from different faces of a cell, it is no problem that the surface field is larger than the volume field; for the flat plate, we only correct normal to the wall; everything else would simply mask the true underlying problem

I think what we have to understand is how to constrain the nuEff values on the first face normal to the wall such that the solution becomes stable. Could you please

1. run a setup at yPus=100 in which only the wall faces are corrected
2. run the case for ~100 iterations
3. write out nuEff for the wall faces and the first faces normal to the wall for every iteration
4. upload the data to the cloud folder

Thanks! Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I uploaded the simulations' data in the cloud folder. There are two cases as follows:

1. Folder "turbulentFlatPlate_wmSF_onlywall_100iters" : Only the wall is corrected.
2. Folder "turbulentFlatPlate_wmSF_1stface_100iters" : The first face is also corrected.

In each "yplus_100" folder, the related data including nuEff is saved as csv files for each iteration step. Therefore, please check the folder in the cloud.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, I had a look at the data, and I believe to have found the issue. First, I noticed that the correction based on Spalding would suggest an increased effective viscosity for the first cell face normal to the wall, which is counterintuitive because for yPlus=100, the normal derivate of the velocity profile should be close to zero. Accordingly, the effective viscosity should be smaller than the molecular one. Therefore, I checked the computation again and found a logical mistake:

You must compute uTau only once based on the cell-centered velocity; this uTau is used for computing both correction terms at the wall as well as at the first face in normal direction; instead, you computed uTau again based on the velocity at the first cell face; moreover, the correct distance to use for the uTau computation is dist_fface[faceI], but this should have no effect in the present setup, since the distance to the face is twice the distance to the center

Moreover, I suggest the following additional output:

• uTau
• the velocity derivatives based on uTau
• the effective viscosity normalized by the molecular one

I know that some quantities could be determined during post-processing, but I would like to have as many consistency checks as possible directly available.

It might be also beneficial to lower the relaxation coefficient for the SIMPLE U update by setting

relaxationFactors
{
    equations
    {
        U               0.3;
        ...
    }
}

in fvSolution.

I hope that these suggestions get us again another step closer to a working implementation :+1: Best, Andre

[JihooKang-KOR]

Dear Andre,

I tested and uploaded two simulations in the cloud folder. The tests are as follows:

1. Correcting until 1st cell face, relaxation factor U 0.3, max iteration steps 500 (turbulentFlatPlate_wmSF_1stface_500iters_U0.3min): ```nuEff[oppFaceIDs[faceI]] = min(nu, nu_*abs(dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL));``` is applied to the first face, dUdy for wall and face is written out.
2. Correcting only the wall patch, relaxation factor U 0.3, max iteration steps 500 (turbulentFlatPlate_wmSF_onlywall_500iters_U0.3_min): dUdy for wall and face is written out. To print the value dUdy_f_spalding, dUdy_f_spalding[faceI] = -nuEff[oppFaceIDs[faceI]]*magGradUf[faceI]/nu_; is added to the code.

Here are the graphs of Cf and residuals. As seen in the figures, they didn't capture the behavior of the original function.

Meanwhile, I saw the values of dUdy_spalding_face in the csv file and found that they are much less than dUdy_spalding_wall except the very front and end of the plate, whereas dUdy for the first face is only slightly smaller than that for the wall. The calculated uTau values are around 4.0, and thus it seems to be reasonable that the slope at the wall differs from the slope at the first face according to the given Spalding's law graph. (In case of y+ = 100, first cell center height = 0.00052226311, first cell face height = 0.0010445262)

According to the above, I couldn't find any wrong calculation, but the result was unfortunately not good. dUdy_spalding_face for case 2 (only wall correction) is rather very huge (normalized nuEff is bigger than 100), but it converges well.

That is what I found so far, and I hope you can find the uploaded csv files in the cloud to investigate precisely.

Best regards, Jihoo Kang

[AndreWeiner]

Dear Jihoo, thanks for the update. The velocity profile you plotted provides very valuable information. As you can see, the profile's velocity value at y=0.001, where we want to evaluate it, is already higher than the inflow velocity (~69m/s). This is a problem of Spalding's function because the velocity profile approaches infinity as y approaches infinity. The uTau we get based on the integral fit is too high to get meaningful velocity derivatives at the first cell face. I estimated that for the first cells along the plate the yPlus value is actually > 300, which is close to or even past the range we can use it.

Here are two experiments I would like you to run:

1. correct both face values, but set the effective viscosity at the first cell face everywhere to 2*nu_molecular
2. correct both face values, but set the effective viscosity at the first cell face everywhere to 0.5*nu_molecular

I would like to see the general trend of these corrections. Best, Andre

[JihooKang-KOR]

Dear Andre,

I ran two experiments and obtained the results as follows:

1. *0.5nu_molecular** (turbulentFlatPlate_wmSF_1stface_500iters_U0.3_0.5numol) As seen in the figures, Cf still fluctuates and the pressure residual doesn't converge.

2. *2nu_molecular** (turbulentFlatPlate_wmSF_1stface_500iters_U0.3_2numol) Cf cannot still capture the actual behavior, but the pressure residual becomes more stable.

According to the above results, I concluded that nuEff at the first face should be much larger to capture the actual behavior. Therefore, I tried several larger values, and finally found the proper value for y+ = 100 as follows (Here, relaxation factor for U is again 0.9, whereas 0.3 for the above cases):

3. *40nu_molecular** (turbulentFlatPlate_wmSF_1stface_U0.9_40numol)
   • y+ = 100

• y+ = 50

The behavior of y+ =50 doesn't match very precisely, but it is quite ok. Thus, the value of around 40 times of the molecular viscosity at the first cell face can make the simulation stable and the result reasonable (but depending on y+).

The whole simulation folders are already in the cloud folder. Thus, I hope you can see the related files there.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo,

very interesting results you got there. It appears that I was wrong the entire time regarding the expected effect at the first cell face. We are actually evaluating the profile at a yPlus where the linear approximation still significantly underestimates the slope at the first cell face. Moreover, the slope of Spalding's function actually never becomes zero but rather constant as the distance to the wall becomes larger. Therefore, it is unlikely that the predicted effective viscosity becomes smaller than the molecular viscosity. I guess the essential question to answer is how to limit the effective viscosity at the beginning of the plate (or in general for high yPlus values > 200 of the first cell face).

Here are some more experiments:

• correction at both faces without limitation of the effective viscosity
• correction at both faces with nuEff=min(nuEff, 100*nu_molecular) for the first face
• correction at both faces with nuEff = (1-w) nu_molecular+w nuEff, where the weight w is defined as w=sigmoid(0.05(200-yPlusFace)); the sigmoid function is defined as sigmoid(x) = 1/(1+exp(-x)); yPlusFace is the yPlus value evaluated at the first cell face yPlusFace = yFaceuTau/nu_molecular

Best, Andre

[JihooKang-KOR]

Dear Andre,

The results of the simulations are as follows:

1. Correction at both faces without limitation of the effective viscosity (turbulentFlatPlate_wmSF_1stface_U0.9_nolimit)

2. *Correction at both faces with nuEff=min(nuEff, 100nu_molecular) for the first face** (turbulentFlatPlate_wmSF_1stface_U0.9_min100nu)

3. Correction at both faces with nuEff = (1-w) nu_molecular+w nuEff (turbulentFlatPlate_wmSF_1stface_U0.9_sigmoid)

All the cases couldn't capture the behavior of Cf and didn't converge. For the case without limitation, only the velocity converged, but nothing converged for the other cases.

I uploaded them in the cloud folder as I did before. Thank you for reading this comment.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for running the experiments! I made the following observations:

• the slope at the first face is almost equal to that at the wall, at least for the first iterations
• the correction factor for the first face is typically significantly larger than at the wall
• you can see within the first 10 iterations if the simulation diverges
• already at the second iteration, there is a large scatter in the derivative at the first face computed with Spalding's function (violet plusses), which is not present in the other fields (uTau, Up, magGradU, du_dy_wall - see figure below)

Maybe, changing the eps value for computing the finite difference can mitigate this scatter (making epsilon larger). If only the wall faces are corrected, you can see the scatter in the derivate at the first face first increasing and then decreasing:

Moreover, by comparing the velocity profiles obtained with a constant nuEff=40 at the first face, you see how the velocity along the wall stops decreasing slowly after several iterations:

In all other simulations, the velocity keeps dropping, and at some point, the difference in the velocity along the wall is so steep that the oscillations set on; once the oscillations start, it is hardly possible to obtain a converged solution:

I assume that we have to take smaller steps in the iterative process. Otherwise, the velocity towards the end of the plate drops too quickly, which in turn leads to a growing difference in effective viscosity and in the velocity itself.

Here is what I suggest:

• investigate the scatter issue, e.g, by changing the epsilon value (lowering it might reduce the scatter)
• for the sigmoid-limiting, investigate the influence of the constants in the exponent (lowering the large value (200) will switch off the model at lower yPlus values; lowering the prefactor smears out the transition)

Maybe you got more insights from the last experiments based on my thoughts.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I would like to urgently inform of a mistake in the code.

For calculating dUdy_f_spalding[faceI], I used Newton's method to find the velocity for all region. However, the value diverges when we use Newton's method for a large cell height. At the first cell face, y+ value is around 300, and therefore I should have written "BISECTION" in the code.

I changed the code as follows and will upload here soon.

forAll (patchFaceIDs, faceI)
{
    if (dist_fface[faceI] < 5e-4)
    {
        dUdy_f_spalding[faceI] = (spalding_velocity(dist_fface[faceI] - eps, uTau[faceI], nu_, "NEWTON", kappa)
        - spalding_velocity(dist_fface[faceI] + eps, uTau[faceI], nu_, "NEWTON", kappa))/(2*eps);
    }

    else
    {
        dUdy_f_spalding[faceI] = (spalding_velocity(dist_fface[faceI] - eps, uTau[faceI], nu_, "BISECTION", kappa)
        - spalding_velocity(dist_fface[faceI] + eps, uTau[faceI], nu_, "BISECTION", kappa))/(2*eps);
    }    
}

As I did in the calculation of uTau, I added if statement in terms of the distance of the first face. Afterward, I could find that dUdy_f_spalding[faceI] values were quite small as you expected in the first place. I'm still not sure if it is now correct or not, but in accordance with the below figure, it might be reasonable that the slope at the first face should be much smaller than that at the wall. (In case of y+ = 100, first cell center height = 0.00052226311, first cell face height = 0.0010445262)

In addition, please see the uploaded surfaceFields csv file here. face_100yp_surfaceFields_iter_284.csv

After changing the code, the values of derivative are less than 10000 for almost everywhere (compared to the wall around 1200000). Some derivative values splash at some points, but the trend is generally less than 10000 (Or greater than -10000).

However, there will be a problem because your previous thinking is correct as per the corrected code. According to this, the normalized nuEff_face is too small. Thus, when we use the min-max function that you provided in the meeting, it will work as the case only with wall correction. And we cannot use the original correction because it cannot capture the behavior.

Thank you for reading this urgent news and I will carefully check what I implemented next time.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, great that you found this bug! Here are two blending functions that should have some effect:

// option 1 - take the maximum of both values
//    nuEff[oppFaceIDs[faceI]] = max(nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL), nuEff[oppFaceIDs[faceI]]);
// option 2 - blend between original and model value; start switching model off beyond yPlusFace 200
    scalar yPlusFace = dist_fface[faceI]*uTau[faceI]/nu_;
    scalar w = 1/(1 + Foam::exp(-0.05*(200 - yPlusFace))); // Sigmoid function
//    nuEff[oppFaceIDs[faceI]] = w*nuEff[oppFaceIDs[faceI]] + (1-w) * nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

Both functions converged for the yPlus=100 case. There are still some odd values in the slope for the first few cells, but they are not as extreme as before. They might be related to the fact that the velocity profile is slightly overshooting in the second cell normal to the wall at the beginning of the plate. I hope we can ignore it for now.

Edit: I suggest also writing out the values for nut, yPlusFace, and the weight w if it is used.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I uploaded a jupyter notebook for wmSimpleFoam model comparision (commit 77e6054). According to the plots in the notebook, we can conclude as follows:

1. Option 1 (max function) cannot be used for very small y+ (= 0.05 and 1). Particularly for y+ = 0.05, there is absolutely no graph.
2. Option 2 (sigmoid function) is the best option so far, but it is not significantly different from the only wall corrected model except y+ = 100.

I think using sigmoid function is reasonable because the concept might be as follows:

• For small y+, large w is used (very small influence of spalding slope)
• For large y+, small w is used (huge influence of spalding slope)

It can be considered as a smooth transition between two values (previous nuEff and current spalding slope). However, if we use another function for w (e.g. parabolic, abs(-logx)+c or any quadratic functions), we are able to get better results, but I'm afraid if this function loses the concept and meaning of using the spalding slope correction.

I uploaded three test cases in the cloud as well, but I wrote out only the last iteration because if I write out fully, the size will be about 3.0GB for one case. Please find the cases turbulentFlatPlate_wmSF_onlywall, turbulentFlatPlate_wmSF_1stface_sigmoid_option1, and turbulentFlatPlate_wmSF_1stface_sigmoid_option2.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks again for the update and the notebook. The effect is not yet as I hoped it would be. The current blending is stable, but it actually switches off the model off where it should be used, namely yplus < 200. I've tried another even simpler approach this morning:

nuEff[oppFaceIDs[faceI]] = 0.5*nuEff[oppFaceIDs[faceI]] + 0.5 * nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

I came to this idea because I realized that nuEff without correction is actually converging towards the value resulting from the wall function when using a blending between the two values. I also think that the pure averaging between the two values only stabilizes the solution algorithm but has a minor influence on the final result. E.g., I would expect a blending of 0.8()+0.2() to deliver similar results.

I ran quick tests for yPlus=30/50/100 and found the following trends for uTau and the correction factor at the first face:

As you can see, the results look better though not yet ideal. There are also some strange wiggles in the curves whose origin I could not yet identify. I have to ponder a little more about the remaining difference. One potential reason might be the limiting of uTau values in the iterative process, which affects typically the first few cells. The maximum uTau value reached by the iterative algorithm depends on the mesh resolution and varies with the mesh resolution normal to the wall.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I tested two simulations as follows:

1. w = 0.5

   nuEff[oppFaceIDs[faceI]] = 0.5*nuEff[oppFaceIDs[faceI]] + 0.5 * nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

2. w = 0.8

   nuEff[oppFaceIDs[faceI]] = 0.8*nuEff[oppFaceIDs[faceI]] + 0.2 * nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

According to the results, case 1 yielded the better graphs than case 2 for greater than y+ = 2, whereas there is no difference for the other y+ values (please refer to the commit 5bc4648). However, case 1 couldn't achieve the better performance than the original wall function model. Even the best behavior exists at y+ = 30, 50, and 100, but they are almost the same with the wall function version.

Therefore, I tried another w coefficient for y+ = 50 and 100 as follows:

3. w = 0.35

nuEff[oppFaceIDs[faceI]] = 0.35*nuEff[oppFaceIDs[faceI]] + 0.65 * nu_*(-dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

As seen in the above figure, the behavior is better than the other 2 cases. I have executed this simulation for all the y+ values in the cluster, and the result will be given later if there aren't any mistakes. I hope it works well.

I uploaded case 1 and case 2 in the cloud as well. Please refer to the folders if you need any information in the csv files.

<Update on 26.08.2021> I made the commit f944d12 to show the result of the case w = 0.35. Almost all the results are better than the other cases. However, it couldn't give the proper result at y+ = 0.05 (also the same for w = 0.5) because magGradU at the face is too small when the mesh is too resolved. Therefore, w might have to be greater than 0.8 (0.8 - 1.0) when y+ is less than 1. Other than this range, we can use w = 0.35.

For y+ = 5 and 10, we still have something to improve, the giggling in the graph. I guess that it is caused by uTau. Both of the graphs are smooth.

The graphs seem to have the similar shape.

As seen the above, uTau and Cf are closely related. Thus, I think the Cf graph will be improved if we can soften the fluctuation by employing any methods.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear @AndreWeiner,

I revised the calculation part of nuEff and tested, but the simulation unfortunately didn't work well. It couldn't converge like the previous simulations for which no modification was applied. Therefore, I set two limitations as follows:

• nuEff at the face is corrected after 50th iteration step in order to stabilize the velocity field.
• Use under relaxation for nuEff by referring to the previous nuEff value.

In addition, I changed the uTau and dUdy_f_spalding calculation part to use only Newton's method due to the following reasons:

• When the initial value of the velocity is changed to magUp_facei = 100*ut, Newton's method also shows the good behavior as Bisection method beyond the first face region at y+ = 100 (compared to IntegralFittingMethod.ipynb, the difference can be seen).

• Yielded dUdy_f_spalding (at face) values by Newton's method are more stable than those by Bisection method at least for several y+ values. When uTau decreases, dUdy_f_spalding should always decrease as well. However, I could find that dUdy_f_spalding values sometimes fluctuated by Bisection method.

I made the commit 09df035 for the above revision. Hence, please refer to the commit if you would like to check something.

The result plots are in the commit 2409827, but you won't find that the results are remarkably different from the case of only wall correction and even more unstable. I guess the reasons for it are as follows:

1. When we see the csv files for each y+, dUdy_sp_face values are written out. For some y+, dUdy_sp_face gradually decreases as uTau decreases. However, for y+ = 1, 5, 10, and 50, there are some unstable values (e.g. uTau = 3.08383125792 -> dUdy_sp_face = -15612.4277565, but uTau = 3.08237916243 -> suddenly dUdy_sp_face = -20357.2323454). This phenomenon influences the downstream and leads to the instability which means the calculation of the slope is not perfectly stable.

2. The ratio of the slopes dUdy_f_spalding[faceI]/magGradUf[faceI] is always greater than 1 for the higher y+ values (greater than the reference curve), whereas around 1 for the smaller y+ (almost the same with the reference curve). When we don't use any relaxation factor for nuEff, the simulation is easy to diverge for the higher y+ since nuEff at the face is quite large (strong turbulence viscosity) which leads to very small magGradUf (at face) due to the very small velocity gap between the first cell center and the face. In case of using any relaxation factor, the simulation converges, but the result is not much different from the only wall correction case because the effect of face correction is reduced by that factor.

I uploaded those two test cases in the cloud (turbulentFlatPlate_wmSF_count50_relaxation0.3_onlyNewton and turbulentFlatPlate_wmSF_count50_relaxation0.5_onlyNewton). Please kindly refer to the csv files to see if the above guesses are reasonable or not.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear @AndreWeiner,

I uploaded the test result in the cloud folder. (turbulentFlatPlate_wmSF_count50_100iter)

In the test folder, there are 3 cases for y+ = 30, 50, and 100, respectively. The correction at the first face starts from 51st iteration step, and each iteration until 100 is saved in each y+ folder.

In addition, I added one more column called Slope_ratio which corresponds to abs(dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL) in the csv files. Then, we are able to see how the slope ratio changes. According to the result, the slope ratio does not change significantly before the correction at the face. However, after 51st iteration, the ratio at the middle and the back of the plate starts to increase. But I need to figure out why it happens.

Thank you for reading this comment, and please refer to the csv files that I uploaded in the cloud.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, I had a look at the data. I would say that the results do not look too bad. After switching on the model, the velocity profile converges within the next 50 iterations (compare the difference between eg. the 90th and 100th iteration). The overall tolerance of 10^-5 might not be reached because there are still small oscillations due to oscillations in the slope ratio. If we mitigate the oscillations in the slope ratio further, the final residual should also decrease. For yplus=50, there are two kinks in all quantities coming from the analytical profile. The kinks are likely related to the iterative process for finding uTau.

For 0.05<x<0.2, the velocity and all derived quantities (also utau, slopes, etc.) are somehow plateauing. This effect is more pronounced as yplus goes up. I also noticed that the velocity normal to the wall at the first cell center is negative in this particular region. I remember that these velocity oscillations were also present without using any modified code (at the very beginning of the plate, the velocity magnitude in the second cell normal to the wall is actually larger than the free stream velocity). As the mesh is refined, these oscillations become smaller. You can quite easily see the oscillation in paraview. I think varying the divergence scheme is not helping the situation if I remember your tests correctly. Therefore I suggest moving the domain boundary to x=0 such that the domain starts where the plate starts. The shorter domain should constrain the velocity at the beginning of the plate a little bit more and might mitigate the oscillations.

As we discussed in the last meeting, additional normalization of utau in the iterative process might stabilize the slope ratio. Moreover, computing the finite difference slope with three points should also yield smoother results, but I would first try the other two steps.

If I didn't express myself clearly enough, let me know and I elaborate more on the ideas. As I mentioned at the beginning, the overall model behavior points in the right direction :-)

Best, Andre

[JihooKang-KOR]

Dear Andre,

I changed the starting point from x = -1 (only in blockMesh, in real it is -0.333...) to x = 0 in the blockMesh, and did the same simulation again which is already uploaded in the cloud (turbulentFlatPlate_wmSF_count50_100iter_x0). In addition, I uploaded the notebook that contains graphs of slope ratio, Ux, Uy, and dU/dy spalding at the first face (with the commit 9cd1408).

According to the notebook, the following points could be found:

1. Regarding the plateau at 0.05 < x < 0.2, it doesn't disappear although I changed the starting point to x = 0. Only difference is that Ux is a bit greater than the original case (x from -1).

2. Regarding the minus value of Uy, it still exists as well. I think the plateau and the minus Uy values appear in the same region. They cannot be mitigated by changing the mesh starting point, and therefore another mitigation idea should be come up with.

3. When we see the graph of the dU/dy spalding slope, we can conclude that the overflow of the slope ratio at y+ = 50 is due to the discontinuity at some points, whereas the other two cases have the smooth curves. Therefore, we can use the case of y+ = 30, and 100 for vanishing the oscillation at the front of the plate before the normalization is implemented.

4. Except Ux at 50 iteration, both cases have the plateau (or negative peak) of slopes, and the negative Uy values. This means that any additional mitigation measures are to be applied regardless of the application of the face correction.

The second thing is about the normalization of the spalding_velocity function to find the velocity from uTau. Unlike a typical normalization of any data, we should change both of a uTau input and a distance input in order to make the exact same shape of spalding's law graph after the normalization. However, I'm not sure if it is possible to get the exact the same velocity value because when we change the scale of uTau and distance within 0 to 1 (or any larger scale), we cannot expect how the shape of the spalding graph changes. Unfortunately, I couldn't come up with any idea so far.

Lastly, I would like to ask two things as follows:

• Regarding the sentence "For yplus=50, there are two kinks in all quantities coming from the analytical profile.", I couldn't find two kinks at y+ = 50. I plotted Cf, but the shape has no kink at y+ = 50 and that of both cases (x = -1 and 0) are similar as follows. Therefore, I would like to see where they are.

• With regards to the sentence "After switching on the model, the velocity profile converges within the next 50 iterations ...", I would like to know if I understand the meaning correctly. I think that the velocity profile seems to converge before 100 iterations (around 90th step), but the simulation is redundantly proceeded more and makes the velocity profiles identical between the first cell center and the first cell face (exceed the limitation of the slope ratio). If it is correct, we should find the way how to stop this redundant calculation by observing the mitigation method of the oscillation.

Thank you for reading this comment about the result.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear @AndreWeiner,

In addition to the above, I tried additional investigations to find where the problem is located.

• As I mentioned before, the simulation goes to overflow due to very small magGradUf (at the face). This is because the difference between the velocity at the 1st cell center and the 1st cell face goes almost 0 (they are almost the same).
• This behavior happens due to very high nuEff at the face compared to nuEff at the cell center (It intrinsically happens because of the face correction). While nuEff at the cell center remains almost the same scale, nuEff at the face goes very high. It means that the viscosity at the face is too high, such that the velocity at the face is limited to the value of the cell center velocity. We can even see that the velocity at the cell center goes faster than the previous time step, whereas the velocity at the face decreases. Consequently, these two velocities become identical.

Thus, I thought that I might be able to modify something at the 1st cell center. Hence, I tried to limit Ux at the cell center which is 95% of Ux at the face as follows:

forAll (patchFaceIDs, faceI)
{    
    // nut already exists at the face, whereas zero at the wall
    nuEff[oppFaceIDs[faceI]] *= abs(dUdy_f_spalding[faceI])/(magGradUf[faceI] + ROOTVSMALL);

    if (U[adjacentCellIDs[faceI]][0] >= 0.95*U_face[oppFaceIDs[faceI]][0])
    {
        U[adjacentCellIDs[faceI]][0] = 0.95*U_face[oppFaceIDs[faceI]][0];
    }
}

There is no more count 50 limitation, and the simulation was done until 100th iteration step for y+ = 30, 50, and 100.

The result was what I expected. The difference between two velocities didn't diminish a lot (due to 5% gap, at least), and magGradUf (at the face) could maintain the reasonable value. The Cf graph for 100th iteration step is as follows:

The valley or kink at Rex 0.2 for y+ = 50 (orange line) is due to the unstable Spalding's slope calculation, which needs to be fixed by normalization. The plateau at the front plate for y+ = 100 (green line) is already mentioned in the above comment, and also should be mitigated. Other than these, the result graph looks quite nice. However, since I artificially limited Ux at the cell center, the residual of Ux couldn't converge (remains the scale of 10^-2) as follows:

Therefore, the following methods might be considered.

1. Keep the method of Ux limitation (95%, or 99%, ...) In this case, the convergence of Ux is difficult to achieve as I mentioned above.
2. Use the same correction method (implemented for the wall and the face) at the 1st cell center I would like to try this method. However, we are modifying the solver itself, and thus I couldn't figure out how to modify nuEff at the cell center. For turbulence->nuEff(), the keyword .internalField() or .internalFieldRef() didn't work unlike editing any turbulence models. I couldn't find any methods to directly edit nuEff of volScalarField in the solver level. I need more time to find it, but I'm not sure if there exists any solution.
3. Use other modification methods at the 1st cell center If No.2 doesn't work, we can think about other ideas.

I uploaded the csv data (where the column of nuEff at the 1st cell center is added) of the case 95% Ux for y+ = 30, 50, and 100 in the cloud. Please refer to the data if you would like to check anything more.

<09.09.2021 Update> I found the way to refer to nuEff at the cell center as follows:

volScalarField& nuEff_vol = turbulence->nuEff().ref();

It is referred to as a reference, such that the value can be edited in the solver. However, since the momentum equation only uses nuEff at the surface, any modification at the cell center doesn't influence the simulation. Even when I employed the term fvc::div(turbulence->nuEff() * dev2(T(fvc::grad(U)))) again, no influence occurred. Therefore, I would like to conclude that the method No.2 doesn't work.

Best regards, Jihoo Kang