Code Monkey home page Code Monkey logo

data_driven_wall_modeling's Introduction

Design and implementation of a data-driven wall function for the velocity in RANS simulations

Introduction

A wall function is used to reduce computational cost by employing an analytical profile. However, this approach still yields mesh-dependent results for coarser meshes. The mesh-dependency comes from the difference between the analytical profile and the numerical profile from a simulation as shown in the video.

spalding_wallunit_pres.mp4

The red-dotted lines at the wall and at the first cell face are the actual velocity gradient (slope), and the orange line is the numerical velocity gradient. It is obvious that the slopes are different for coarser meshes. In addition, the actual velocity value (blue graph) and the numerical velocity value (orange graph) at the first cell face are also different for coarser meshes. The wall function approach somehow corrects this discrepancy at the wall, but it does not deal with the difference at the first cell face. Therefore, one of the objectives of this project is to mitigate the discrepancies at the wall as well as the first cell face by correcting the diffusive and the convective fluxes in the momentum equation. Another objective is to create ML models instead of the analytical profile to predict actual values.

Dependencies

  • SingularityCE 3.8.4
  • OpenFOAM-v2106
  • PyTorch 1.9.0 (only CPU)
  • Jupyter-lab 3.2.1
  • Numpy 1.21.3
  • Pandas 1.3.4
  • Matplotlib 3.1.2
  • Flowtorch 1.0

Getting started

Compiling the data-driven wall function and running the test cases requires an installation of OpenFOAM-v2106. Other releases might work as well but have not been tested.

To build the image at a local system, a docker container should be downloaded as follows.

git clone https://github.com/AndreWeiner/of_pytorch_docker.git
cd of_pytorch_docker
docker build -t user_name/of_pytorch:of2106-py1.9.0-cpu -f Dockerfile .

Afterward, this docker image can be converted to a singularity image as follows.

sudo singularity build of2106-py1.9.0-cpu.sif Singularity.def

After the creation of the singularity image, it should be moved to the top folder of this project. If the name of the image is not of2106.sif, please revise the name to of2106.sif, and then the below instructions can be executed.

The more detailed instruction for this version is available in this link written by Dr.-Ing. Andre Weiner. If the latest version is needed, please refer to this repository.

Data generation phase

1D channel

  1. The case is found in ./test_cases/channelFlow1D.
  2. To run the test case, create a run folder (ignored by version control) only for the first time, copy the case from test_cases to run, and execute the Allrun_1D script.
mkdir -p run
cp -r test_cases/channelFlow1D run/
cd run/channelFlow1D
./Allrun_1D
  1. The simulation consists of two parts, 50 time steps with 0.0002 write interval and 90 time steps with 0.001 write interval, but they are automatically executed by the Allrun_1D script.

  2. After the simulations are finished, the 1D data extraction utility is to be used, which is found in (topFolder)/utilities/extractData.

  3. In the folder, type wmake to compile the extraction utility.

  4. Then, move to ../../run/channelFlow1D folder, and type extractData. This makes a csv file that contains all the features and the labels.

  5. In (topFolder)/notebooks, BestMapping_1DChannel.ipynb needs to be executed to obtain the trained ML models. In the data loading cell in this notebook, the folder name needs to be changed if a specific folder name is used in the run folder. In the case of this document, channelFlow1D_1.388_69.4 should be changed to channelFlow1D.

field1d_path = run + 'channelFlow1D_1.388_69.4/extractData.csv'
field1d_data = pd.read_csv(field1d_path, delim_whitespace=False)
  1. After the notebook is executed, the two folders are created to use in OpenFOAM. These models are already in each 2D case if needed.

Application of trained ML models

2D flat plate (No wall function / Standard wall function cases)

  1. The case is found in ./test_cases/turbulentFlatPlate.
  2. To run the test case, copy the case from test_cases to run, and execute the Allrun.singularity script.
cp -r test_cases/turbulentFlatPlate run/
cd run/turbulentFlatPlate
./Allrun.singularity
  1. To change the type of wall function, the nut file in (topFolder)/0.SpalartAllmaras folder needs to be modified. For no wall function, keep the script as it is. For the standard wall function nutUSpaldingWallFunction, comment out the line 45 and 46, and comment the line 48 and 49.
42    bottomWall
43    {
44        // No wall function
45        type            fixedValue;
46        value           uniform 0.0;
47        // Standard wall function
48        //type            nutUSpaldingWallFunction;
49        //value           $internalField;
50    }
  1. This simulation is executed as parallel, and therefore MPI is used. To change the setting of the number of CPU cores, edit decomposeParDict in (topFolder)/system folder.

2D flat plate (Data-driven wall function)

  1. The case is found in ./test_cases/turbulentFlatPlate_datadriven.
  2. To run the test case, copy the case from test_cases to run, and execute the Allrun.singularity_ddwmSimpleFoam_serial script.
cp -r test_cases/turbulentFlatPlate_datadriven run/
cd run/turbulentFlatPlate_datadriven
./Allrun.singularity_ddwmSimpleFoam_serial
  1. To change the correction method (either the wall correction or the wall/face correction), move to (topFolder)/solvers/datadriven_wmSimpleFoam and edit nuEffCorrection.H. For the wall correction, the line 108 and 109 should be commented out, then no face correction occurs. For the wall/face correction, keep it as it is.
105 // Blended 1st face correction with yPlus blending
106 forAll (oppFaceIDs, faceI)
107 {
108     scalar blendFace = nuEff[oppFaceIDs[faceI]]*(labelAccessor[faceI][1]*U_ref/l_ref)/(face_b[faceI] + ROOTVSMALL);
109     nuEff[oppFaceIDs[faceI]] = w_face[faceI]*nuEff[oppFaceIDs[faceI]] + (1 - w_face[faceI])*blendFace;
110 }
  1. wmake needs not to be executed because it is already included in the Allrun.singularity_ddwmSimpleFoam_serial script.

  2. In (topFolder)/notebooks, PlotCf_2DflatPlate.ipynb shows the mesh-dependency of four scenarios that correspond to no wall function, standard wall function, data-driven wall function with wall correction, and data-driven wall function with wall/face correction by plotting skin friction, which was investigated in the thesis. In order to execute the notebook, change the folder name in the notebook accordingly for one's cases.

  3. The folders normalized1DModel and scaleModule made from the ML training phase are already included in the test case, and thus one can run the data-driven wall function cases without executing the notebook BestMapping_1DChannel.ipynb in the training phase.

  4. This simulation is executed as serial due to an indexing issue along the plate in the modified solver.

2D NACA-0012 airfoil (No wall function / Standard wall function cases)

  1. The case is found in ./test_cases/airFoil2D_Re3e6_alpha0. The case is with Re = 3e6 and the angle of attack 0 degree.
  2. To run the test case, copy the case from test_cases to run, and execute the Allrun_serial.singularity script.
cp -r test_cases/airFoil2D_Re3e6_alpha0 run/
cd run/airFoil2D_Re3e6_alpha0
./Allrun_serial.singularity
  1. To change the type of wall function, the nut file in (topFolder)/0.orig folder needs to be modified. For no wall function, keep the script as it is. For the standard wall function nutUSpaldingWallFunction, comment out the line 38 and 39, and comment the line 40 and 41.
36    airfoil
37    {
38        type            fixedValue; // without wall functions
39        value           uniform 0;
40        //type            nutUSpaldingWallFunction; // with the wall function
41        //value           $internalField;
42    }
  1. The NACA profile for blockMesh is from the repository https://github.com/AndreWeiner/naca0012_shock_buffet.

2D NACA-0012 airfoil (Data-driven wall function)

  1. The case is found in ./test_cases/airFoil2D_Re3e6_alpha0_datadriven.
  2. To run the test case, copy the case from test_cases to run, and execute the Allrun_ddwmSFairfoil_serial.singularity script.
cp -r test_cases/airFoil2D_Re3e6_alpha0_datadriven run/
cd run/airFoil2D_Re3e6_alpha0_datadriven
./Allrun_ddwmSFairfoil_serial.singularity
  1. To change the correction method (either the wall correction or the wall/face correction), move to (topFolder)/solvers/ddwmSimpleFoam_airfoil and edit nuEffCorrection.H (for diffusive flux) and fluxCorrection.H (for convective flux). For the wall correction, the line 138 and 139 for nuEffCorrection.H as well as the line 14 and 15 for fluxCorrection.H should be commented out, then no face correction occurs. For the wall/face correction, keep it as it is.

nuEffCorrection.H

133 // Blended 1st face correction with yPlus blending
134 forAll (oppFaceIDs, faceI)
135 {
136     if (mesh.Cf().boundaryField()[surfaceID][faceI].x() < 0.998)
137     {
138         scalar blendFace = nuEff[oppFaceIDs[faceI]]*(labelAccessor[faceI][1]*U_ref/l_ref)/(face_b[faceI] + ROOTVSMALL);
139         nuEff[oppFaceIDs[faceI]] = w_face[faceI]*nuEff[oppFaceIDs[faceI]] + (1 - w_face[faceI])*blendFace;
140     }
141 }

fluxCorrection.H

 9 // Convective flux correction at the 1st face
10 forAll (oppFaceIDs, faceI)
11 {
12     if (mesh.Cf().boundaryField()[surfaceID][faceI].x() < 0.998)
13     {
14         scalar blendPhi = phi[oppFaceIDs[faceI]]*(Uface_b[faceI])/(mag(U_face[oppFaceIDs[faceI]]) + ROOTVSMALL);
15         phi[oppFaceIDs[faceI]] = w_face[faceI]*phi[oppFaceIDs[faceI]] + (1 - w_face[faceI])*blendPhi;
16     }
17 }
  1. In (topFolder)/notebooks, PlotCf_Cp_2Dairfoil.ipynb shows the mesh-dependency of four scenarios by plotting skin friction and pressure coefficient, which was investigated in the thesis. In order to execute the notebook, change the folder name in the notebook accordingly for one's cases.

  2. The folders normalized1DModel and scaleModule made from the ML training phase are already included in the test case as well.

Notebooks

Used in the thesis

  • BestMapping_1DChannel.ipynb : Creation of the ML models that used in the thesis.
  • Mapping_1DChannel.ipynb : Checking uncertainties of each hyper-parameter to find the best setting of hyper-parameters.
  • MeshEffect_Spalding.ipynb : Showing the effect of coarser meshes based on Spalding's function (the notebook for explaining mesh-dependency).
  • PlotCf_2DflatPlate.ipynb : Plot skin friction for the flat plate case at Re = 1e7.
  • PlotCf_2DflatPlate_Re3e6.ipynb : Plot skin friction for the flat plate case at Re = 3e6.
  • PlotCf_2DflatPlate_Re6e6.ipynb : Plot skin friction for the flat plate case at Re = 6e6.
  • PlotCf_Cp_2Dairfoil.ipynb : Plot skin friction and pressure coefficient for the airfoil case at Re = 3e6.
  • PlotCf_bestBlendingExponent_2DflatPlate.ipynb : Finding proper blending exponents gammaWall and gammaFace for diffusive fluxes.
  • FieldAnalysis_avgNuEff.ipynb : Field analysis for average effective viscosity over iteration steps and along the plate for y+ = 0.05, 10, and 100.

Others

  • BestMapping_distanceBasedModel_1D.ipynb : Creation of the ML models for distance based wall model that is similar to the standard wall function.
  • FieldAnalysis_1stFace.ipynb : Field analysis over iteration steps such as nuEff and face slopes to find the reason of divergence before the reverse blending is employed for the 1st cell face correction.
  • PlotCf_distanceBasedModel_2DflatPlate.ipynb : Plot skin friction for the flat plate case at Re = 1e7 with respect to the distance based wall model.
  • PlotResiduals_2Dairfoil.ipynb : Plot residuals to initially check convergence status for the airfoil case at Re = 3e6.
  • PlotResiduals_2DflatPlate.ipynb : Plot residuals to initially check convergence status for the flat plate case at Re = 1e7.
  • Plots_monotonicVelLimit.ipynb : Plot for monotonic velocity limitation setting before the reverse blending is employed for the 1st cell face correction.

Solvers

Used in the thesis

  • datadriven_wmSimpleFoam : The modified solver from simpleFoam for the flat plate cases by employing the trained ML models.
  • ddwmSimpleFoam_airfoil : The modified solver from simpleFoam for the airfoil case by employing the trained ML models.

Other

  • ddwmSimpleFoam_wallCorr : The distance-based modified solver for the flat plate cases.

Utilities

Used in the thesis

  • extractData : Data extraction utility for the 1D channel flow.

Others

  • extractData_2D : Data extraction utility for the 2D flat plate case to compare 1st face slopes between the 1D channel and the 2D flat plate with resolved mesh.
  • extractData_cellCenterVel : Data extraction utility for the 1D channel flow with respect to the distance-based model.

Thesis

The thesis for the project : DOI

BibTeX :

@misc{kang_jihoo_2022_6590747,
  author       = {Kang, Jihoo},
  title        = {{Design and implementation of a data-driven wall 
                   function for the velocity in RANS simulations}},
  month        = may,
  year         = 2022,
  publisher    = {Zenodo},
  doi          = {10.5281/zenodo.6590747},
  url          = {https://doi.org/10.5281/zenodo.6590747}
}

data_driven_wall_modeling's People

Contributors

jihookang-kor avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar

data_driven_wall_modeling's Issues

The cause of skin friction discrepancy for airfoil cases

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.

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

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.

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

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

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

The flux correction at the 1st face

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

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

Slopes for the finest mesh in a 2D flat plate case

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

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.