ProjectPhysX / FluidX3D

The fastest and most memory efficient lattice Boltzmann CFD software, running on all GPUs via OpenCL. Free for non-commercial use.
https://youtube.com/@ProjectPhysX
Other
3.81k stars 301 forks source link

Issue with setting up inflow/outflow boundaries. The boundaries keep reflecting back. #202

Closed He11Bee closed 1 week ago

He11Bee commented 2 months ago

src.zip For a school project i wish to measure the forces of the wind on different versions of the Ahmed body. The simulation works but I am not able to set the inflow/outflow boundaries. The wind keeps bouncing back at the boundaries. See figure. Bouncing [EQUILIBRIUM_BOUNDARIES] is uncommented and the boundaries are set to TYPE_E. if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n] = TYPE_E; Any ideas of what can the problem be.

There is also an issue on the front wall as it creates a strange formation, 1-wall strange formations

This is the code I'm using:

void main_setup() { // Ahmed body; required extensions in defines.hpp: FP16C, FORCE_FIELD, EQUILIBRIUM_BOUNDARIES, SUBGRID, optionally INTERACTIVE_GRAPHICS // ################################################################## define simulation box size, viscosity and volume force ################################################################### const uint memory = 6000u; const float lbm_u = 0.05f; const float box_scale = 2.0f; const float si_u = 60.00; const float si_nu=1.48E-5f, si_rho=1.225f; const float si_width = 0.976021f, si_height = 0.848059f, si_length = 2.61945f; // Geometry for ahmed_25deg_m_Spoiler0 // const float si_width= 1.45359f, si_height= 2.05124f, si_length= 2.61945f; // Geometry for ahmed_25deg_m_Spoiler12 // const float si_width=1.67423f, si_height=1.14702f, si_length=3.69274f; // Geometry for Hatcback_NoSpoiler const float si_A = si_widthsi_height+2.0f0.05f0.03f; const float si_T = 0.25f; const float si_Lx = units.x(box_scalesi_width); const float si_Ly = units.x(5/8box_scalesi_length); const float si_Lz = units.x(box_scalesi_height); const uint3 lbm_N = resolution(float3(1.67423f, 3.69274f, 1.14702f), memory); units.set_m_kg_s((float)lbm_N.y, lbm_u, 1.0f, box_scalesi_length, si_u, si_rho); print_info("Re = "+to_string(to_uint(units.si_Re(si_width, si_u, si_nu)))); const float lbm_nu = units.nu(si_nu); const float lbm_length = units.x(si_length); LBM lbm(lbm_N, lbm_nu); // ###################################################################################### define geometry ###################################################################################### const string model = "ahmed_25deg_m_Spoiler12"; // input the name of the model, so that i can create new .txt files ever time i run a new sim to save the data. Mesh* mesh = read_stl(get_exe_path()+"../stl/" + model + ".stl", lbm.size(), lbm.center(), float3x3(float3(0, 0, 0), radians(0.0f)), lbm_length); mesh->translate(float3(0.0f, 0.0f, 2.0f - mesh->pmin.z)); lbm.voxelize_mesh_on_device(mesh, TYPE_S|TYPE_X); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z); if (z==0u) lbm.flags[n] = TYPE_S; if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = lbm_u; if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n]=TYPE_E; }); // ####################################################################### run simulation, export images and data ########################################################################## lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_Q_CRITERION; lbm.graphics.set_camera_centered(0.0f, 0.0f, 80.0f, 1.648722f); lbm.run(0u); // initialize simulation

if defined(FP16S)

// const string path = get_exe_path()+"FP16S/"+to_string(memory)+"MB/";

elif defined(FP16C)

const string path = get_exe_path()+"FP16C/"+to_string(memory)+"MB/";

else // FP32

// const string path = get_exe_path()+"FP32/"+to_string(memory)+"MB/";

endif // FP32

//lbm.write_status(path);
//write_file(path+"Cd.dat", "# t\tCd\n");
const string path2 = "C:/..../Exports/" + model + ".txt";
while(lbm.get_t()<=units.t(si_T)) { // main simulation loop
    for (int i = 0; i < 10; i++) // "i < number", where the number determined the number of times that I collect data.
    {
        lbm.run(100u); // run lbm_dt LBM time steps
        lbm.calculate_force_on_boundaries(); 
        lbm.F.read_from_device();
        lbm.voxelize_mesh_on_device(mesh, TYPE_S);
        const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S | TYPE_X);
        print_info("Step nr : " + to_string(lbm.get_t()));
        print_info("SI - Drag: " + to_string(units.si_F(lbm_force.y)) + " N " + "SI - Lift: " + to_string(units.si_F(lbm_force.z)) + " N");
        write_line(path2, "Step nr : " + to_string(lbm.get_t()) + "\t" + "SI - Drag: " + to_string(units.si_F(lbm_force.y)) + " N" + "\t" + "SI - Lift: " + to_string(units.si_F(lbm_force.z)) + " N" + "\n");

if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)

//      lbm.graphics.write_frame(path+"images/");

endif // GRAPHICS && !INTERACTIVE_GRAPHICS

    }
    lbm.run(1u);
}   
//lbm.write_status(path);

} /**/

cTatu commented 1 week ago

I replicates the same scenario. I downloaded the latest release version 2.18 and just run the ahmed body example from source.cpp as is. GIF of fluid reflecting from boundary walls (which were supposedly set to TYPE_E): ezgif com-webp-to-gif-converter

Is this desired default behavior? Also the Drag coefficient doesn't seem to converge even after many steps: Program output:

|                                     \ /               FluidX3D Version 2.18 |
|                                      '     Copyright (c) Dr. Moritz Lehmann |
|-----------------------------------------------------------------------------|
| Info: Unit Conversion: 1 cell = 6.650 mm, 1 s = 180460 time steps           |
| Info: Re = 1577027                                                          |
|----------------.------------------------------------------------------------|
| Device ID    0 | NVIDIA H100                                                |
| Device ID    1 | NVIDIA H100                                                |
| Device ID    2 | NVIDIA H100                                                |
| Device ID    3 | NVIDIA H100                                                |
|----------------'------------------------------------------------------------|
| Info: OpenCL C code successfully compiled.                                  |
| Info: Allocating memory. This may take a few seconds.                       |
| Info: Loading                                                               |
|       "/gpfs/scratch/bsc19/bsc19959/FluidX3D/bin/../stl/ahmed_25deg_m.stl"  |
|       with 28626 triangles.                                                 |
| Warning: Some boundary cells have non-zero velocity, but MOVING_BOUNDARIES  |
|          is not enabled. If you intend to use moving boundaries, uncomment  |
|          "#define MOVING_BOUNDARIES" in defines.hpp.                        |
|-----------------.-----------------------------------------------------------|
| Grid Resolution |                              702 x 1884 x 189 = 249965352 |
| Grid Domains    |                                             2 x 2 x 1 = 4 |
| LBM Type        |                                    D3Q19 SRT (FP32/FP16C) |
| Memory Usage    |                               CPU 6913 MB, GPU 4x 4040 MB |
| Max Alloc Size  |                                                   2264 MB |
| Time Steps      |                                                     45115 |
| Kin. Viscosity  |                                                0.00000185 |
| Relaxation Time |                                                0.50000556 |
| Reynolds Number |                                             Re < 58833072 |
|---------.-------'-----.-----------.-------------------.---------------------|
-38.530 0.300
-13.483 0.285
-2.020 0.287
1.642 0.322
0.305 0.447
3.284 0.514
-1.012 0.400
3.862 0.331
-0.251 0.321
2.960 0.335
0.103 0.455
2.022 0.336
1.562 0.425
1.558 0.379
1.084 0.360
-0.061 0.288
0.785 0.301
-0.105 0.365
1.207 0.356
-0.257 0.425
1.443 0.414
-0.204 0.343
1.532 0.311
0.146 0.416
1.268 0.327
0.868 0.480
1.033 0.298
1.086 0.344
0.207 0.320
1.039 0.343
0.251 0.309
1.025 0.330
0.301 0.317
1.148 0.296
0.201 0.285
1.161 0.330
0.789 0.302
1.111 0.430
0.630 0.294
0.948 0.291
0.829 0.344
0.353 0.303
0.991 0.415
0.282 0.315
0.951 0.280
0.277 0.298
1.105 0.364
0.332 0.301
1.017 0.347
0.675 0.475
0.991 0.312
0.711 0.310
0.938 0.312
0.756 0.297
0.399 0.342
0.863 0.302
0.336 0.311
0.951 0.312
0.289 0.326
1.045 0.322
0.322 0.372
0.963 0.291
0.633 0.368
0.892 0.344
0.740 0.358
0.866 0.478
0.768 0.317
0.430 0.334
0.802 0.322
0.434 0.320
0.885 0.333
0.382 0.339
0.989 0.310
0.347 0.303
0.961 0.380
0.577 0.330
0.891 0.319
0.674 0.325
0.833 0.340
0.757 0.484
0.438 0.280
0.790 0.299
0.438 0.302
0.818 0.317
0.412 0.301
0.943 0.341
0.356 0.294
0.935 0.302
0.578 0.320
0.914 0.352
0.614 0.332
0.842 0.316
0.701 0.325
0.443 0.327
0.769 0.306
0.401 0.339
1.004 0.475
0.403 0.312
0.964 0.465
0.380 0.314
0.918 0.316
0.575 0.295
0.904 0.314
0.596 0.316
0.875 0.318
0.656 0.299
0.453 0.292
0.760 0.303
0.386 0.317
1.003 0.325
0.365 0.312
1.009 0.315
0.384 0.454
0.938 0.323
0.598 0.333
0.882 0.299
0.661 0.311
0.851 0.327
0.686 0.337
ProjectPhysX commented 1 week ago

Hi @He11Bee @cTatu,

it is expected that pressure waves weakly reflect off TYPE_E boundaries. The TYPE_E cells themselves "eat"/discard all incoming DDFs and replace them with equilibrium DDFs with specified density/velocity. But the layer of regular fluid cells next to a TYPE_E wall gets from one side incoming DDFs from a shockwave, and on the other the equilibrium DDFs from the TYPE_E cells, and they collide. In this adjacent regular fluid layer, some compression can occur, leading to shockwave reflection.

Forces for the Ahmed body will never converge - it is a transient solution with the vortices explicitly resolved, so forces will always fluctuate. The time average will converge though.

Kind regards, Moritz

cTatu commented 1 week ago

Thank you very much for the explanation!

He11Bee commented 1 week ago

Thank you for the replication and answer.


From: Dr. Moritz Lehmann @.> Sent: Tuesday, September 3, 2024 7:11:33 AM To: ProjectPhysX/FluidX3D @.> Cc: He11Bee @.>; Mention @.> Subject: Re: [ProjectPhysX/FluidX3D] Issue with setting up inflow/outflow boundaries. The boundaries keep reflecting back. (Issue #202)

Hi @He11Beehttps://github.com/He11Bee @cTatuhttps://github.com/cTatu,

it is expected that pressure waves weakly reflect off TYPE_E boundaries. The TYPE_E cells themselves "eat"/discard all incoming DDFs and replace them with equilibrium DDFs with specified density/velocity. But the layer of regular fluid cells next to a TYPE_E wall gets from one side incoming DDFs from a shockwave, and on the other the equilibrium DDFs from the TYPE_E cells, and they collide. In this adjacent regular fluid layer, some compression can occur, leading to shockwave reflection.

Forces for the Ahmed body will never converge - it is a transient solution with the vortices explicitly resolved, so forces will always fluctuate. The time average will converge though.

Kind regards, Moritz

— Reply to this email directly, view it on GitHubhttps://github.com/ProjectPhysX/FluidX3D/issues/202#issuecomment-2325614131, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BJXLNJ6APBUADDSZZVWFVFTZUVAILAVCNFSM6AAAAABKQ6QWLSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMRVGYYTIMJTGE. You are receiving this because you were mentioned.Message ID: @.***>