jihookang-kor / data_driven_wall_modeling Goto Github PK

Dear @AndreWeiner,

In addition to the 1st issue, I would like to mention about the correction at the 1st face.

In the 1st issue, I showed that the 2D mapping model didn't work well because the slopes at the wall were too large to converge. Meanwhile, the slopes at the final step (3539th) for the finest mesh case were the similar to the range of the slopes for the 1D channel flow case. Therefore, I kept the original 1D mapping model for the simulation, and came up with another approach.

Regarding the notebook "MeshEffect_Spalding.ipynb" (the commit 997d829), the graph in the notebook can be applied only when the velocity profile of the real simulation is similar to the Spalding's function profile. However, when we turn on the face correction, the velocity profile drastically changes, and thus it isn't similar to the Spalding's function anymore. In contrast, the wall correction works well (steadily increasing wall shear stresses), and smoothly changes the velocity field and the velocity gradient at the wall.

Hence, I decided that I fix "magGradUf" (velocity gradient at the 1st face) at a certain time step as the commit 8d809a7. In the code, I fixed the magGradUf at 6th time step with the variable name of "fixedGradUf". I was not sure if the velocity profile at this time step was similar to the Spalding's function, but I thought that this was the time step before the drastic change in the velocity gradient at the face. Therefore, I tried this setting for y+ = 30, 50, and 100, and the result of skin friction is as follows:

The result shows that the simulation at least converged. If we let the face correction happen with changing magGradUf every time step, the magGradUf diminishes very fast, and then the correction ratio of nuEff always exceeds 1. This makes the simulation crashes because magGradUf reduces further and further.

However, when the face correction with the fixed denominator (fixedGradUf) at a certain early time step turns on, the nuEff correction ratio is less than 1 which can lead to convergence. Therefore, my speculation is that we can properly correct fluxes at the 1st face as well if we are able to find which time step has the velocity profile that is similar to the Spalding's function.

I need to think about any methods how to find it if my speculation is reasonable.

Best regards,
Jihoo Kang

Dear @AndreWeiner,

With reference to the subject, I made a code to extract wall slopes and face slopes in a 2D flat plate case (commit 1bab3b1). Then, I extracted the related data from the case y+ = 0.05. Subsequently, I made two ML models for a correction test as follows.

1. Wall slope model
   a. features : x-coordinates, Ux at 1st face
   b. label : wall slope

2. Face slope model
   a. features : x-coordinates, y-coordinates at faces, integral average Ux
   b. label : face slope

At first glance, the slope values in the 2D case look similar to the values in the 1D channel case when the 2D case converged (at 3539th step). However, at the beginning of the simulation, the slopes at the wall are much larger, whereas the slopes at the faces are much smaller than the range of the 1D case. Therefore, when the 2D ML models were applied to the correction, an over-correction happened at the wall, an under-correction happened at the 1st face, compared to using the 1D mapping data.

When we use the 1D mapping data or the slope calculation of Spalding's function, the correction at the 1st face strongly exerts on the face, and thus the gap of the velocities between the 1st cell center and the 1st cell face is almost zero (then, the simulation crashes).

The situation for the 2D ML model is, however, opposite. In this case, the correction at the wall is too strong, so the simulation explodes at the wall.

I guess the slopes in the 1D mapping only fit for near the last time step, and hence the model cannot reflect the correct slopes at the beginning of the simulation. Therefore, I think the amount of the correction should gradually be decreased for the wall, and be increased for the 1st face. For instance, 2.0*wallslope at the beginning, then changes gradually to 1.0*wallslope after the certain steps, while 0.2*faceslope at the beginning, then to 1.0*faceslope. But it is only a rough and vague idea.

Please refer to the files and folders in the "2DMapping" folder in the cloud as follows.

1. The notebook "2DFlatPlateMapping.ipynb"
   Contains two 2D mapping models in the notebook.
2. The folder "turbulentFlatPlate_noWallFunc_ref0.05_datasets"
   Slope datasets with the finest mesh simulation in the folder. "extractData_wall_2D.csv" and "extractData_face_2D.csv" are the dataset csv files.
3. The folder "turbulentFlatPlate_2Dmodel"
   The simulation by using the 2D mapping model.

Thank you very much for reading this issue.

Best regards,
Jihoo Kang

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