Trouble with force readouts

[Maere05]

Hi,

First of all, thank you for this amazing software. The performance is out of this world. I'm trying to simulate a wingsuit at different angles of attack. I was using Ansys Fluent for this but it just takes forever on a PC at home and a limited license, even at a way lower resolution.

I'm having trouble with the force readouts, expected lift and drag at an AoA of 10° is in the hundreds of Newtons. (That’s why I commented out the scientific notation, since 500 is easier to read than 5.0E2). I get 55N drag and 15N lift if I don’t convert to SI, or 6.62E7 N drag and 1.74E7 N lift with units.si_F();

I assume that I made mistake in the setup of the conversion between LBM and SI units. And I find weird that the drag is way bigger that the lift, wingsuits should have a glide ratio of around 3.

Do you see any obvious mistakes in the setup.cpp ?

Cheers and thank you for your work! Marius

Setup.cpp

void main_setup() {

    const uint res = 240;               // Grid Resolution
    const float AoA = 10.0f;            // Angle of Attack [°]

    // -- SI units --
    const float si_Length = 2.0f;       // Characteristic length[m]
    const float si_Airspeed = 50.0f;    // Airspeed in SI [m/s]
    const float si_rho = 1.225f;        // 1.225kg/m^3 air density
    const float si_nu = 1.48E-5f;       // 1.48E-5f m^2/s kinematic shear viscosity of air

    // -- LBM units --
    const float lbm_x = 1.0f;
    const float lbm_u = 0.1f;
    const float lbm_rho = 1.0f;

    units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_Length, si_Airspeed, si_rho);

    LBM lbm(res*1.2, res*1.6, res*1.2, units.nu(si_nu));

    const float size = 1.0f * (float)res;
    const float3 center = float3(lbm.center().x, lbm.center().y, lbm.center().z); // center the model
    const float3x3 rotation = float3x3(float3(1, 0, 0), radians(30.0f-AoA)) * float3x3(float3(0, 0, 1), radians(00.0f)) * float3x3(float3(1, 0, 0), radians(90.0f)); // Stl is Y-up -> Rotate by 90° | 30.0f because the model is already tilted

    lbm.voxelize_stl(get_exe_path() + "stl/wingsuit.stl", center, rotation, size, TYPE_S); // Import Mesh

    const uint N = lbm.get_N(), Nx = lbm.get_Nx(), Ny = lbm.get_Ny(), Nz = lbm.get_Nz();

    for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < N; n++, lbm.coordinates(n, x, y, z)) {
        if (lbm.flags[n] != TYPE_S) {   
            lbm.u.y[n] = lbm_u;     //All non Solid nodes get the Velocity
        }
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) {
            lbm.flags[n] = TYPE_E;  // Flag Inlets and Outlets
        }
    }

    key_1 = true;// Show Mesh
    key_4 = true;// Show Turbulence

    lbm.run(1500);// Run    

    for (int i = 0; i < 10; i++)    // Run for 100 Steps, then print the Forces to the Console
    {
        lbm.run(100);

        lbm.calculate_force_on_boundaries();
        lbm.F.read_from_device();
        const float3 force = lbm.calculate_force_on_object(TYPE_S);
        print_info("Step nr : " + to_string(lbm.get_t()));
        print_info(/*"Lateral: " + to_string(force.x) +*/ "LBM - Drag: " + to_string(force.y) + " N Lift:" + to_string(force.z) + " N");
        print_info(/*"Lateral: " + to_string(force.x) +*/ "SI  - Drag: " + to_string(units.si_F(force.y)) + " N Lift:" + to_string(units.si_F(force.z)) + " N");
    }

    wait(); // wait for a keypress to close the Program
}

Console Out:

` |----------------.------------------------------------------------------------| | Device ID 0 | NVIDIA GeForce GTX 1080 | |----------------'------------------------------------------------------------| |----------------.------------------------------------------------------------| | Device ID | 0 | | Device Name | NVIDIA GeForce GTX 1080 | | Device Vendor | NVIDIA Corporation | | Device Driver | 522.25 | | OpenCL Version | OpenCL C 1.2 | | Compute Units | 20 at 1797 MHz (2560 cores, 9.201 TFLOPs/s) | | Memory, Cache | 8191 MB, 960 KB global / 48 KB local | | Buffer Limits | 2047 MB global, 64 KB constant | |----------------'------------------------------------------------------------| 1 warning generated. | Info: OpenCL C code successfully compiled. | Loading "C:/dev/FluidX/bin/stl/wingsuit.stl" with 203691 triangles. | Info: Voxelizing mesh. This may take a few minutes. | |-----------------.-----------------------------------------------------------| | Grid Resolution | 288 x 384 x 288 = 31850496 | | LBM Type | D3Q27 SRT (FP32/FP32) | | Memory Usage | CPU 516 MB, GPU 3796 MB | | Max Alloc Size | 3280 MB | | Time Steps | 1500 | | Kin. Viscosity | 0.00000001 | | Relaxation Time | 0.50000004 | | Reynolds Number | Re < 2644988928 | | Volume Force | 0.00000000, 0.00000000, 0.00000000 | |---------.-------'-----.-----------.-------------------.---------------------| | MLUPs | Bandwidth | Steps/s | Current Step | Time Remaining | | Info: Step nr : 1600 | | Info: LBM - Drag: 62.980694 N Lift:11.393983 N | | Info: SI - Drag: 77151344.000000 N Lift:13957628.000000 N | | Info: Step nr : 1700 | | Info: LBM - Drag: 61.000553 N Lift:11.799897 N | | Info: SI - Drag: 74725672.000000 N Lift:14454872.000000 N | | Info: Step nr : 1800 | | Info: LBM - Drag: 59.884335 N Lift:13.046824 N | | Info: SI - Drag: 73358304.000000 N Lift:15982358.000000 N | | Info: Step nr : 1900 | | Info: LBM - Drag: 59.182091 N Lift:14.444882 N | | Info: SI - Drag: 72498056.000000 N Lift:17694978.000000 N | | Info: Step nr : 2000 | | Info: LBM - Drag: 58.040073 N Lift:15.479063 N | | Info: SI - Drag: 71099080.000000 N Lift:18961850.000000 N | | Info: Step nr : 2100 | | Info: LBM - Drag: 57.142467 N Lift:15.655202 N | | Info: SI - Drag: 69999520.000000 N Lift:19177620.000000 N | | Info: Step nr : 2200 | | Info: LBM - Drag: 55.727894 N Lift:15.647128 N | | Info: SI - Drag: 68266664.000000 N Lift:19167730.000000 N | | Info: Step nr : 2300 | | Info: LBM - Drag: 55.218250 N Lift:15.055589 N | | Info: SI - Drag: 67642352.000000 N Lift:18443094.000000 N | | Info: Step nr : 2400 | | Info: LBM - Drag: 53.899487 N Lift:14.339896 N | | Info: SI - Drag: 66026864.000000 N Lift:17566372.000000 N | | Info: Step nr : 2500 | | Info: LBM - Drag: 54.105560 N Lift:14.259664 N | | Info: SI - Drag: 66279304.000000 N Lift:17468086.000000 N | | 1082 | 235 GB/s | 34 | 2500 100% | 0s |

`

[trparry]

For an external flow problem, you might want to create a larger fluid domain. Fluid domain size is an important source of error for external flow where boundary conditions close to the solid object could cause problems. Just an idea based on when I compare theoretical results to simulation for external flow over a cylinder.

[RedBlaze42]

Have you achieved good results with the cylinder ?

In that case, can you share your main_setup ?

[trparry]

@RedBlaze42 It's probably easier to adapt the "cylinder in rectangular duct" setup already in setup.cpp, but I used a 0.0254 m dia x 0.25m length cylinder from an stl file because I wanted to observe the flow effects around the ends.

Theoretical drag force was somewhere between 0.388N and 0.46N in the flow direction (y for me), depending on the drag coefficient you take (drag coeff should be around 1.0-1.2 based on the reynold's number). Simulation reported 0.57N.

The delta could be a number of different factors including grid resolution, cylinder length (theoretical calc is for an infinitely long cylinder), the duration of the simulation (I only did 2s), or perhaps I still need to increase the size of the fluid domain relative to the diameter to reduce the effects of the BC's. My grid resolution was actually quite low for the cylinder size that I used (the stair step effect on the face of the curved cylinder is probably throwing off my results). I should move to a higher grid resolution, but its a decent first validation.

You can reproduce the stl easily, the aspect ratio between length and diameter is the only thing that matters.

void main_setup() {
    // make a list of all variables you have in SI units (m, kg, s)
    const float si_x = 0.25f; // 0.25m cylinder length
    const float si_u = 10.0f; // 10m/s wind speed
    const float si_rho = 1.225f; // 1.225kg/m^3 air density
    const float si_nu = 1.48E-5f; // 1.48E-5f m^2/s kinematic shear viscosity of air

    // set velocity in LBM units, should be between 0.01-0.2
    const float lbm_u = 0.1f;
    const float lbm_rho = 1.0f; // density in LBM units always has to be 1
    const float lbm_x = 256u; 

    // set unit conversion factors between SI units and LBM units
    units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_x, si_u, si_rho);

    // compute kinematic shear viscosity in LBM units
    const float lbm_nu = units.nu(si_nu);

    // set grid resolution based on lbm_x
    const uint Nx = to_uint(1.25*lbm_x);
    const uint Ny = to_uint(lbm_x);
    const uint Nz = to_uint(lbm_x/2);

    LBM lbm(Nx, Ny, Nz, lbm_nu); // create LBM object

    // load geometry from stl file, mark all grid points that belong to cylinder geometry with (TYPE_S|TYPE_X)

    const float size = 1.0f * (float)lbm_x;
    float3 center = float3(lbm.center().x, lbm.center().y, lbm.center().z);
    const float3x3 rotation = float3x3(float3(0, 0, 1), radians(90.0f));
    Mesh* mesh = read_stl("<path to your stl file>", lbm.size(), center, rotation, size);
    mesh->center = float3(lbm.center().x, 0.25f * size, lbm.center().z); 
    voxelize_mesh_hull(lbm, mesh, TYPE_S);

    // set box boundary conditions
    for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < lbm.get_N(); n++, lbm.coordinates(n, x, y, z)) {
        //if (!(lbm.flags[n] & TYPE_S)) lbm.u.y[n] = lbm_u; // initial velocity
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E;
        //if (z == 0u) lbm.flags[n] = TYPE_S;
        if (lbm.flags[n] != TYPE_S) lbm.u.y[n] = lbm_u;
    }

    // run simulation
    lbm.run(units.t(2.0f)); // run for 2 seconds

    // calculate forces on boundaries on GPU, then copy force field to CPU memory
    lbm.calculate_force_on_boundaries();

    lbm.F.read_from_device();

    // sum force over all boundary nodes marked with (TYPE_S)
    const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S);

    // transition back to SI units
    const float si_force_x = units.si_F(lbm_force.x); // force in Newton
    const float si_force_y = units.si_F(lbm_force.y);
    const float si_force_z = units.si_F(lbm_force.z);

    // print x/y/z force components
    print_info("z force = " + to_string(si_force_z) + " N");
    print_info("y force = " + to_string(si_force_y) + " N");
    print_info("x force = " + to_string(si_force_x) + " N");
}/**/

[trparry]

Also, you should look at the Stokes drag setup in setup.cpp for another validation with theory.

[RedBlaze42]

Sorry for being a little out of subject, I can create another issue if you want.

I'm having trouble with reproducing your results. I'm using a sphere because I can find easily the Cd on wikipedia (0.47 for a Reynolds of 10^4)

I'm not sure about the length I use in the reynolds calculation. And also I don't know the impact of the lbm speed on the results because it should be canceled by the unit conversion.

Have you an idea of what I'm missing ?

void main_setup() {
    // make a list of all variables you have in SI units (m, kg, s)
    const float si_x = 0.25f; // 0.25m sphere radius
    const float si_u = 0.2975f; // 0.2975m sphere radius to adjust for Reynolds number
    const float si_rho = 1.225f; // 1.225kg/m^3 air density
    const float si_nu = 1.48E-5f; // 1.48E-5f m^2/s kinematic shear viscosity of air

    // set velocity in LBM units, should be between 0.01-0.2
    const float lbm_u = 0.05f; // what should I use ?
    const float lbm_rho = 1.0f; // density in LBM units always has to be 1
    const float lbm_x = 480u;
    const float R = 150.0f;

    // set unit conversion factors between SI units and LBM units
    units.set_m_kg_s(R, lbm_u, lbm_rho, si_x, si_u, si_rho);

    // compute kinematic shear viscosity in LBM units
    const float lbm_nu = units.nu(si_nu);

    print_info("Re= "+to_string(units.si_Re(2.0f*si_x, si_u, si_nu))); // Needs to be 10^4 to match wikipedia's 0.47 Cd

    // set grid resolution based on lbm_x
    const uint Nx = to_uint(lbm_x);
    const uint Ny = to_uint(lbm_x);
    const uint Nz = to_uint(lbm_x);

    LBM lbm(Nx, Ny, Nz, lbm_nu); // create LBM object

    const float size = 1.0f * (float) lbm_x;
    float3 center = float3(lbm.center().x, lbm.center().y, lbm.center().z);

    // set box boundary conditions
    for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < lbm.get_N(); n++, lbm.coordinates(n, x, y, z)) {
        //if (!(lbm.flags[n] & TYPE_S)) lbm.u.y[n] = lbm_u; // initial velocity
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E;
        //if (z == 0u) lbm.flags[n] = TYPE_S;
        if(sphere(x, y, z, lbm.center(), R)){
            lbm.flags[n] = TYPE_S;
        }else lbm.u.y[n] = lbm_u;
    }

    // run simulation
    lbm.run(5000u); // run for 5000 timesteps
    key_4 = true;
    lbm.graphics.set_camera_centered((float) (lbm.get_t()/37), 25.0f, 100.0f, 0.95f);
    lbm.graphics.write_frame_png(get_exe_path()+"export/rotation/");
    lbm.graphics.set_camera_centered(0.0f, 90.0f, 90.0f, 0.5f);
    lbm.graphics.write_frame_png(get_exe_path()+"export/t/");

    // calculate forces on boundaries on GPU, then copy force field to CPU memory
    lbm.calculate_force_on_boundaries();

    lbm.F.read_from_device();

    // sum force over all boundary nodes marked with (TYPE_S)
    const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S);

    // transition back to SI units
    const float si_force_x = units.si_F(lbm_force.x); // force in Newton
    const float si_force_y = units.si_F(lbm_force.y);
    const float si_force_z = units.si_F(lbm_force.z);

    // print x/y/z force components
    print_info("z force = " + to_string(si_force_z) + " N");
    print_info("y force = " + to_string(si_force_y) + " N");
    print_info("x force = " + to_string(si_force_x) + " N");
}```

[Maere05]

Ill try with the bigger bounding box, but i think then the will be a problem. Maybe a bigger GPU will be the solution... Is there another way to "Fix" this? because my values are of by a factor of 10, the Cylinder is of by 4%. Thanks!

[RedBlaze42]

I don't think that it's only the bounding box becaus @trparry achieved great results with a 256*320*128 box.

[trparry]

@RedBlaze42, you are using a low velocity. See the setup from the author in issue #16 where he uses si_u = 3.0f; // 3m/s wind speed yet lbm_u = 0.1f; also see #27 for his comments on units.

Let's assume a constant grid resolution of 256x320x128, my code then scales the solid object within this grid. A large object that fills most of this grid will have more voxels that describe it's curved surface, which is good. However, this large object will also be closer to the boundaries of the domain, which is bad. On the other hand, if I scale the object to be very small, only a few voxels would describe the surface of the object (bad), but it would be far away from the boundaries (good). So, its a tradeoff for a limited amount of computing power. If you have a more capable GPU, you could just make a very large grid of 1 billion voxels and have both an object that is far away from the boundaries, but yet is described by a sufficient number of voxels. Indeed, this is exactly what some of the larger test cases are doing, they use grids that are >1billion voxels.

I observed different results if my object was close to the boundaries. However, as I said I think grid resolution could be a contributing factor as well. I would try increasing your grid resolution, and scaling down your object and compare the results to when the object is larger/closer to the boundaries. @Maere05 I would start with a very simple test case with a simple geometry before jumping to a complex geometry first to make sure everything works.

[trparry]

@RedBlaze42 I also noticed that the forces oscillate wildly between positive and negative at the start of the simulation, where you need to run for a sufficient number of timesteps in order to wait for the forces to stabilize. This behavior makes sense based on how you set the velocity in your flow field. To visualize how these force vectors are oscillating in windows, toggle the 1 key.

[RedBlaze42]

@trparry For my tests I now run 10k timesteps. I achieved good results for the Cx of a sphere and a cube when they're small compared to the simulation box. But I still don't understand what length I should use for the Reynolds calculation. And what lbm_speed is ideal. Is it arbitrary ?

Thanks for all your help

[RedBlaze42]

As I said I got great results for basic volumes with 0.1 lbm_speed and a very large sim box compared to the volume. Here I have the graph of the coefficient of drag of a thin disc according to the timesteps and it's never stable even at 50k+ timesteps. Is that expected ? (sorry I'm new to CFD)

Here is the code https://pastebin.com/HTviZu7j I'm using FP16C and SRT

[Maere05]

@trparry Thanks, Ill do some testing on a simple airfoil shape, the sphere drag checks out. how big was your boundary in comparison to the object? In Ansys Fluent I always used +100% on all 6 sides of the object bounding box. Ill do some testing in FluidX3D, maybe LBM needs bigger bounds. But at some point the resolution of the object will just get too small since I only have 8gb of VRAM right now. Cheers

[trparry]

@Maere05 you may need to do a sensitivity analysis on a simple test case, where you keep the size of a voxel relative to the object the same but vary the size of your domain (will involve math), at some point the flow field won't change much as you increase the domain size. In my experience rules of thumb overestimate the size of the domain needed to reasonably approximate an unbounded domain.

@RedBlaze42 just to clarify is this a plot of the drag coefficient? There aren't axis labels. I would think your plot is reasonable for drag coeff, since this is presumably highly turbulent flow (I don't know your re number). Turbulence is by definition an unsteady phenomena, as vortices detach from the object this will create local areas of pressure fluctuation on the surface, resulting in a noisy drag force which is what we would expect in the physical world. The length used for reynolds number is the characteristic length, there are many resources online to learn about how to choose this. Concerning the velocity in lbm units, lbm uses a different unit system from SI units, to have the numerical values close to 1.0f where accuracy and precision are best. So, I think choosing 0.1 for velocity enables better accuracy and precision.

[RedBlaze42]

@trparry Sorry for not being clear, yes the y axis is drag coefficient and x is timesteps. With the length of the cube as the characteristic length I have a 5E4 Reynolds number.

Thank you for explaining the lbm units. For the characteristic length my teacher said that I sould use the total length of the object along the airspeed axis when calculating the drag coeff is that correct ?

I'm making a script to test different box to object ratio do you want me to share the results here ?

[Maere05]

Ok, so I did lots of Simulations of a NACA 2412 airfoil in a 40 m/s airflow. The wing is always 64 voxels wide, and the sim runs for 50k timesteps. Then I change the size of the bounds. I started at 96 and worked up to a boundary that is 6x bigger than the wing at its widest. These results are very confusing, has anyone had similar experiences?

[trparry]

@Maere05 an object that is only 64 voxels wide might not be able to represent a curved airfoil accurately enough to do a meaningful benchmark. Could you show a screenshot of the flagged type s voxels (the 1 key) so we can see it? Also maybe share your setup. If you need a better GPU so that you can do larger simulations, maybe consider using a GPU hosted on Amazon aws or the like.

[Maere05]

Hi, I also thought 64 is not a lot, so I did a new run with 128.

The code I used for the different scales is quite basic :

void main_setup() { // Airfoil

    for (int iter = 1; iter < 12; iter++)
    {
    //int iter = 8;
        float boundaryScale = 1.0f + (0.25f * (float)iter);
        print_info(to_string(boundaryScale));

        const uint L = 128u;
        const float size = 1.0f * (float)L;
        const float si_x = 1.0f; // 1m chord length
        const float si_u = 40.0f; // 40m/s wind speed
        const float si_rho = 1.225f; // 1.225kg/m^3 air density
        const float si_nu = 1.48E-5f; // 1.48E-5f m^2/s kinematic shear viscosity of air

        // set velocity in LBM units, should be between 0.01-0.2
        const float lbm_u = 0.1f;
        const float lbm_rho = 1.0f; // density in LBM units always has to be 1
        const float lbm_x = L;

        // set unit conversion factors between SI units and LBM units
        units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_x, si_u, si_rho);

        // compute kinematic shear viscosity ion LBM units
        const float lbm_nu = units.nu(si_nu);

        // set grid resolution based on lbm_x
        const uint Nx = to_uint(L * (1.0f * boundaryScale));
        const uint Ny = to_uint(L * (2.0f * boundaryScale));
        const uint Nz = to_uint(L * (1.0f * boundaryScale));

        LBM lbm(Nx, Ny, Nz, lbm_nu); // create LBM object

        // #############################################################################################################################################################################################
        const float3 center = float3(lbm.center().x, lbm.center().y - 30, lbm.center().z);
        const float3x3 rotation = float3x3(float3(1, 0, 0), radians(-4.0f)) * float3x3(float3(0, 0, 1), radians(90.0f)) * float3x3(float3(1, 0, 0), radians(90.0f));
        lbm.voxelize_stl(get_exe_path() + "/stl/airfoil-naca-2412.stl", center, rotation, size, TYPE_S);// | TYPE_X)

        const uint N = lbm.get_N();
        for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < N; n++, lbm.coordinates(n, x, y, z))
        {
            if (lbm.flags[n] != TYPE_S) lbm.u.y[n] = lbm_u;
            if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E; // all non periodic
        }
        //double E1 = 10.0, E2 = 10.0, E0 = 10.0;

        print_info("");
        print_info("New sim Starts Here ! ");
        print_info("Nx = " + to_string(Nx));

        for (int i = 0; i < 500; i++)
        {
            lbm.run(100);
            lbm.calculate_force_on_boundaries();
            lbm.F.read_from_device();
            const float3 force = lbm.calculate_force_on_object(TYPE_S);// | TYPE_X
            //print_info("Lateral: " + to_string(units.si_F(force.x)) + "N - Drag: " + to_string(units.si_F(force.y)) + " N - Lift:" + to_string(units.si_F(force.z)) + " N");
            print_info("Step nr:" + to_string(lbm.get_t()) + " Lift: " + to_string(units.si_F(force.z)) + " Drag:" + to_string(units.si_F(force.y)));
        }
    }
}

And here is a Screenshot of the Airfoil

[ProjectPhysX]

Hi @Maere05 @trparry @RedBlaze42, just briefly skimming through and commenting. Thanks for the great discussion here! Getting these validation test cases to work is not easy. From my experience, about 1-3x the size of the object (here length of the airfoil front-to-back) in either direction (top-bottom and front-back) works best for the compromise of under-resolved geometry and boundary effects. For the boundaries left-right of the airfoil, it's better to make them not TYPE_E but periodic (no boundaries) instead, because then the boundaries will not directly touch the object, and the airfoil effectively has infinite length. You can also decrease the width of the box then (for turbulence testing it shouldn't be too small however, maybe 32 nodes at least) so there is more resolution for the other 2 dimensions.

The force oscillations after initialization is due to pressure waves bouncing off the object when the initially moving fluid hits it head-on. After a few tousand time steps (more for larger resolution) they should be gone, absorbed by the non-reflecting TYPE_E boundaries.

[trparry]

@ProjectPhysX thank you! Ok yes the pressure oscillations after initialization makes sense. So, the larger the resolution the longer you may have to wait for the these oscillations to dissipate, because it takes more time steps for these pressure waves to reach the type E boundaries due to the larger domain. @Maere05 I would try to add a "waiting" period before calling lbm.calculate_force_on_boundaries();, might have to run a few tests to figure out how many steps to wait. Perhaps try your case at 160 because this is where the oscillations become more severe, such that you can try to wait for 50k or 100k steps but it won't take forever in real world time.

print_info("");
        print_info("New sim Starts Here ! ");
        print_info("Nx = " + to_string(Nx));
lbm.run(XXXX); //waiting period dependent on resolution
        for (int i = 0; i < 500; i++)
        {
            lbm.run(100);
            lbm.calculate_force_on_boundaries();
            lbm.F.read_from_device();
            const float3 force = lbm.calculate_force_on_object(TYPE_S);// | TYPE_X
            //print_info("Lateral: " + to_string(units.si_F(force.x)) + "N - Drag: " + to_string(units.si_F(force.y)) + " N - Lift:" + to_string(units.si_F(force.z)) + " N");
            print_info("Step nr:" + to_string(lbm.get_t()) + " Lift: " + to_string(units.si_F(force.z)) + " Drag:" + to_string(units.si_F(force.y)));
        }

[RedBlaze42]

I did the opposite of @Maere05's tests: I fixed the box size at 620 and tested different cube/sphere sizes.

I made two plots per object (I didn't forget the legends this time !) The first plot is the drag coefficient by timesteps repeated for every object size. The second plot is the average drag coefficient for every box_size/object_size ratio, the average is taken from timestep 16k to 40k.

The simulations are at 10^4 reynolds so 0.025m sphere radius or 0.05 cube length, 0.2975 m/s wind speed, 1.225 rho and 1.48E-5 m²/s kinetic viscosity. FP32/FP16C and SUBGRID activated because the simulation exploded without subgrid, even at 10^4 Reynolds and I didn't see any difference using FP32/FP32 (good job @ProjectPhysX !!)

For lbm units: lbm_u = 0.125 and lbm_rho=1.0

What do you think about this results ?

[Maere05]

Very interesting results! I did another test of the Airfoil, this time without the periodic boundaries in the X axis. The object was 1/12th of the Z axis. I did it to 1 million Timesteps, but the oscilations did not change. (the Lift makes a cool patter tho :) )

Il continue testing with a lower object/bounds ratio. I hope to find a good balance between resolution - timesteps - accuracy.

@RedBlaze42 have you ever had oscilating forces like the 1/32 that converge later on? In my experience, after 10k steps, when the sim has settled in and all the turbulence is "where it should be", the results don't change much.

[RedBlaze42]

@Maere05 the small objects never converges. The 1/32 ratio means that the cube/sphere is 20 voxels wide, I think that it's too small to have a stable force.

The drag coefficient match the one on wikipedia when the box is 6 times larger than the object.

[RedBlaze42]

For a school project I used FluidX3D and found that when taking the average force over 15k timesteps, with multiple object/box ratio, when there is a plateau on the graph (like the 2nd and 4th graph on my precedent post), the result as an error of 10% compared to Ansys Fluent but we don't know yet if FluidX3D is closer to the real value than Fluent is. Fluent also takes into consideration the roughness of the material which FluidX3D doesn't I think.

[Maere05]

Hi Everyone, I continued to test on the Airfoil section. The force lines you get when pressing 1 were very stable at lower resolutions, and were "jiggly" at higher resolutions, until they were just all over the place. Changing the length of the fluid domain had very little effect, it was pretty much just the number of lattice points between the airfoil and the top of the domain. So using a higher resolution airfoil made it less stable.

https://user-images.githubusercontent.com/62503321/203380964-acdd086d-6f41-47d1-b2de-a410b0605423.mp4

https://user-images.githubusercontent.com/62503321/203380966-2b7d4ea6-db57-435c-a9c5-75e3a48d15a4.mp4

The issue with a very narrow path between the airfoil and the wall is, that it will constrict the airflow more and create a different low pressure zone above the wing. To see if the forces would ever stabilze, i simulated variations from Nz=64 to Nz=512 to 300k steps. after 50k steps I did not see any change.

The last one i continued to 2M steps, but that also didn't not have an effect.

The L/D ratio is still way too low, according to airfoiltools.com it should be around 80. Ill continue testing on higher resolutiuon airfoils, as these are still quite blocky. Is there a better method than exporting the forces to a .txt file? ill try to do a basic UI that plots the forces live, but that will take some time.

[trparry]

@Maere05 you said that "force lines you get when pressing 1 were very stable at lower resolutions, and were "jiggly" at higher resolutions" but looking at your plots the amplitude of the noise on the lift decreases, or at least stays relatively constant as resolution increases. Yes, the drag is more noisy with higher resolution. A few additional questions:

1.) This looks to be a nearly a 2D domain now. How many voxels wide is your domain? It should be at least 32 voxels wide as mentioned, probably should be more. 2.) I assume the sides of your domain are non periodic (no boundaries), and the top, bottom, inlet, and outlet are all type E correct? 3.) Is your airfoil still 128 voxels? Or, since you've moved to 2D why can't you make it much higher resolution? Why stick to Nz=64 to Nz=512? Since now the width of your domain is so much smaller, you can increase the size of the other two dimensions, which also increases the resolution of the airfoil as well.

[Maere05]

What I ment with the jiggly lines is, that at higher resolutions they will suddenly point into the airfoil even when there is a low pressure zone above it (even at the start of the foil, not just where there is flow separation).

Correct: X is non periodic, Z and Y are type E

For these sims, the chord was 128 units long, I saw that they looked more stable at the lower resolutions, so I tested there first. A sim at a chord length of 512 (128X2048X1024) just finished and the results are a similar. Average Drag: 7.7N Average Lift: 97.9N L/D = 12.3

I found that at lower resolutions, the simulations look more stable and the direction of the forces intuitively makes more sense. But the average over +100k steps is more accurate. Ill continue testing and hope to find something close to the 80/1 L/D ratio.

[trparry]

@Maere05 there may be something wrong with your theoretical benchmark. As mentioned here airfoiltools.com is strictly 2D, where induced drag is not considered. Your simulation is not truly 2D, where the number of voxels between the two periodic boundaries is presumably about 32 or more voxels. This means your simulation may include induced drag (which is substantial at low air speeds), while your theoretical benchmark does not. I'm not entirely sure about this though. You'll have to do more research to verify that what you're simulating actually matches your theoretical benchmark.

Also, I found a paper here that used Fluent on NACA 2412 and they got an L/D of 12 for a 7 deg angle of attack and an airspeed of 31.89 m/s. This is very close to your results. I would read the paper and look at how they setup this simulation. You are making the assumption that CL/CD = L/D which I think is correct but maybe there is something more subtle going on here.

Another approach to validation is to use an FVM navier stokes based code like Fluent (costs money) or OpenFoam (free). You could also try Simscale which is also free and most of it uses OpenFoam in the background but has an easier user interface (incidentally, they also have an LBM solver). Validating between LBM and FVM navier stokes is just a good sanity check, and I think many here would be interested in the results (particularly the differences in results between well known standard test cases where many physical experiments have been done and it is widely accepted in the community. Also comparisons between compute time is helpful because most FVM NS solvers in most cases would be slower than FluidX3D). I think some of the test cases already compared navier stokes with LBM, but its nice to see different use cases because LBM has its own niche that it is valid for (most notably for air, Mach < 0.3 because you can make the incompressibility approximation). There are ways to use LBM for compressible flows but this has not been implemented in FluidX3D (yet).

[trparry]

Also, it is worth saying again that voxel-based brick shaped meshing will detrimentally effect the physical accuracy of the simulation due to the stair step effect, compared to meshes with tetrahedral elements that can conform better to arbitrarily curved surfaces. Brick shaped meshes might be faster than tetrahedral with a GPU, or maybe tetrahedral meshes are more difficult to implement with LBM, @ProjectPhysX can probably comment more on this. But, if we increase the resolution of the brick shaped voxel grid high enough, the stair step effect starts to become negligible. The point at which this happens depends on the setup and your requirements (mesh sensitivity study). So, when you compare your results to standard/theoretical test case geometry like an airfoil or cylinder, you have to keep this in mind when speaking about accuracy between simulation and theory. If you are getting tripped up by 1-5% delta between simulation and theory, this very well might be the culprit. If you are wanting to get within sub 1% of theory, you should use a very high resolution grid, and the computation time might be quite long. Also, the physical experiments underlying standardized test cases like drag on a cylinder have a tolerance on repeatability themselves so this must be taken into account as well.

Side note: it would be interesting to use a standard test case (cylinder) and compare tetrahedral meshing + navier stokes FVM vs FluidX3D + brick shaped voxel meshing where the meshes have similar resolutions. Or even play around with grid refinement on a navier stokes FVM code and compare compute time and the difference in results.

[Maere05]

Thank you for your notes on Airfoiltoos, that makes a lot of sense. I did a comparison of Fluent (LES) vs FluidX3D. Both airfoils had a chord length of 1m and are 1m wide at an AoA of 8°. The results were somewhat close. Fluent got me 40N of drag and 275N of Lift (L/D = 6.8). FluidX3D got me 76N of drag and 210N of lift (L/D = 2.7). When increasing the resolution in FluidX3D, the L/D keeps rising and getting closer. Which intuitively makes sense, bricks are not very aerodynamic. I think ill upgrade my GPU, because this is filling my 8gb of VRAM. The other way of getting these sims closer would probably be an adaptive lattice, which @ProjectPhysX mentioned in his masters theisis is planned to be implemented in the future. This is how Fluent Meshing does it: But im sure that will be a pain to implement and is way outside of my own capabilities. Fluent was quite a bit faster, makes sense, since the resolution is vastly different (500k vs 150M). Ill continue testing on a bigger GPU and on higher AoAs since that is what im interested in most.

[RedBlaze42]

I think ill upgrade my GPU, because this is filling my 8gb of VRAM..

I can only suggest using kaggle's free notebooks equipped with a Quadro P100 with 16GB of VRAM. When FluidX3D will be able to use multiple GPUs at once you can even use the notebook with 2x Tesla T4 with both 16GB of VRAM for a total of 32GB. They're a bit slower than the P100 but there's two of them and it's still free. I checked and using kaggle for other things than AI is still in there TOS. The setup needs a bit of thinkering but it may be easier than upgrading your GPU.

[Maere05]

Thanks for all your help to everyone! When I max out the resolution on my new RTX 3090 (so worth it :) ) It matches the values from the paper @trparry mentioned. I get an L/Dof 11 at AoA 8°. So to anyone who has trouble with forces in the future, use higher Resolution. Cheers Marius

[randomwangran]

Side note: it would be interesting to use a standard test case (cylinder) and compare tetrahedral meshing + navier stokes

Check this out: https://github.com/OpenFOAM/OpenFOAM-dev/tree/05ffb6a6ff64a90bb3f0a11da3dfedfee792556c/tutorials/basic/potentialFoam/cylinder