The cause of skin friction discrepancy for airfoil cases

[JihooKang-KOR]

Dear @AndreWeiner,

I would like to report some issues as follows.

1. Interpolation scheme: I tried to write interpolate(U) bounded Gauss linearUpwind grad(U); in fvSchemes, but for the interpolation scheme, I could use linearUpwind or linearUpwindV as similar schemes. When I wrote one of them, the simulation didn't work with the following message.

--> FOAM FATAL IO ERROR: (openfoam-2106 patch=211215)
attempt to read beyond EOF

file: /home/jihookang/0_Projects/Data_driven_wall_modeling/run/airfoil_ddwmSimpleFoam/airFoil2D_Re3e6_alpha0_ddwmSF_upwindTest/system/fvSchemes.interpolationSchemes.interpolate(U) at line 47.

    From virtual Foam::Istream& Foam::ITstream::read(Foam::token&)
    in file db/IOstreams/Tstreams/ITstream.C at line 478.

Therefore, I think we should use linear as the interpolation scheme.

2. Correction of phi: I used the new field phiSgs for the convective flux correction, but the result was the identical to the original one. Hence, I just returned it to the typical field phi. In order to check if the convective flux correction works or not, I multiplied the factor 3.5 to the blending term as,

   scalar blendPhi = 3.5*phi[oppFaceIDs[faceI]]*(Uface_b[faceI])/(mag(U_face[oppFaceIDs[faceI]]) + ROOTVSMALL);

   and the result is as follows:

We can see that it fluctuates, and thus it means that the correction at least works. Afterward, I checked the convective flux correction ratio at the 1st face when I did not apply any additional factor to the blending term as follows.

scalar blendPhi = phi[oppFaceIDs[faceI]]*(Uface_b[faceI])/(mag(U_face[oppFaceIDs[faceI]]) + ROOTVSMALL);

The range of the ratio is within 0.9 ~ 1.15. It means that the convective flux correction in this case doesn't influence much. Instead, the reason why the graph is not matched with the reference one will be explained in No.3.

3. Diffusive flux correction at the 1st face: I once explained that the face correction increased skin friction for flat plate cases. However, the increase happened only at the middle and the end of the plate. At the front plate, the correction factor was very small, and we could find the unstable kink at the front as well. The same thing happens in airfoil cases, at the front and at the end of the airfoil this time as follows:

If the correction factor (ratio) is less than 1, the skin friction values decrease. Particularly, they significantly decrease at the front because the ratio is really small. Therefore, I turned on the face correction only if the correction ratio is larger than 1.0, and the result is as follows:

In this case, only wall correction is applied at the front airfoil. Then, it looks better, but there is still discrepancy. I guess it is the identical phenomenon to the kink in flat plate cases. This uncertainty seems not be able to be solved currently.

The reason of the discrepancy is quite obvious. I would like to decide the method (current blending or no face correction less than 1.0 ratio?). However, this ratio limitation method only works well for the higher y+. I'm afraid if this new method breaks the balance in flat plate cases (especially at y+ = 5 or 10, I additionally would like to add that the face correction ratio for small y+ (less than y+ = 10) for airfoil cases is mostly larger than 1.0 (regardless of the region), and therefore the shape of the graph is totally different from the graphs for the higher y+).

Otherwise, we use the current blending method for all the cases, and we can just report the phenomena of Cf discrepancy caused by the data-driven approach.

Thank you for reading this issue.

Best regards, Jihoo Kang

[AndreWeiner]

Dear Jihoo, thanks for the update and analysis. The correct way of defining the interpolation schemes should be

interpolate(U)  linearUpwindV grad(U);

Gauss is only needed for divergence terms, and bounded adds an additional term to the transport equation, which is also related to the divergence.

Before I can comment on the face blending, there is another issue. What I totally forgot is that the phi depends on the face velocity, too (this is different from the scalar transport case). In the incompressible case, phi is simply the cell-centered velocity vector interpolated to the face and multiplied (dot-product) with the face normal vector. Therefore, phi should be overwritten as follows:

scalar blendPhi = Foam::sign(phi[oppFaceIDs[faceI]]) * Foam::pow(Uface_b[faceI], 2) / (mag(U_face[oppFaceIDs[faceI]]) + ROOTVSMALL);

However, this correction will not work, because, for the phi prediction, we would need the velocity vector and not just the speed. Therefore, I think it is not possible at the moment to implement the convective flux correction properly since this would require restarting from scratch.

Regarding the face correction of the diffusive fluxes, a modification of the yplus-based blending might do the trick. Currently, you are switching it off as yplus becomes relatively small. You can do the same for large values, e.g., like this.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

Thank you for your comment.

I would like to ask two things about your comment as follows:

1. You mentioned that Phi correction using speed would not work, and thus we should restart from scratch. I would like to ask if this means I currently have to implement a new version of phi correction from the beginning.

2. Regarding the diffusive face correction, you proposed a different y+ blending method for airfoil cases. However, I'm afraid if the consistency of the concept could be damaged because the shape of the y+ blending function for airfoil cases will be totally different from the function for flat plate cases. So far I used the function that has only the different coefficient, but if I use the new function, the blending concept will change and it would be hard to find the proper reason why this blending is also needed for higher y+.

I will apply interpolate(U) linearUpwindV grad(U); from now on, but the simulation will be executed after the above things are clear.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo,

1) please do not start from scratch; keep the implementation as it is right now; correcting phi will be a future task (for somebody else) 2) if we change the blending, the blending should be consistent for both setups; my suggestion was to add an upper bound for the y-plus blending factor such that only coarse meshes and the beginning of the plate/airfoil are affected, but maybe underestimate the impact

Best, Andre