Forces explode after Revoxelisation

[Maere05]

Hi everyone, I would like to iterate on an airfoil shape by revoxelising it during a running sim, so that for small changes I don't need as many timesteps to get a result. Restarting it for a tiny change would be a lot slower. To experiment I tested quite large changes (+ 5° AoA) and the airflow looks perfect after some time has passed. I do this by setting everything that is TYPE_S at the moment to be a fluid, then everything that is the new airfoil is TYPE_S and all its values get reset. Then I write everything to the GPU and run more steps.

The airflow looks perfect, like I mentioned, the forces are all over the place. In all the places where it was a fluid at first but is now a solid, it shows a ridiculous amount of drag, so I reduced the scaling factor. I tried everything I thought could help resetting the forces, but I always get the same results. This is how I try to reset the force:

lbm.voxelize_stl(get_exe_path() + "/stl/NACA-Wide.stl", center, rotation, size, TYPE_S);

lbm.run(1000);

        for (int i = 0; i < 10; 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(to_string(lbm.get_t()) + " " + to_string(units.si_F(force.z)) + " " + to_string(units.si_F(force.y)));
        }

        lbm.u.read_from_device();
        lbm.F.read_from_device();
        lbm.rho.read_from_device();
        lbm.flags.read_from_device();

                // remove the old airfoil
        for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < N; n++, lbm.coordinates(n, x, y, z)) // remove the old airfoil
        {
            if (lbm.flags[n] == TYPE_S)
            {
                lbm.u.x[n] = 0.0f;
                lbm.u.y[n] = 0.0f;
                lbm.u.z[n] = 0.0f;
                lbm.rho[n] = 0.0f;

                lbm.flags[n] = TYPE_N;  // 0b00000000
            }
        }
                rotation = float3x3(float3(1, 0, 0), radians(-5.0f)) * float3x3(float3(0, 0, 1), radians(90.0f)) * float3x3(float3(1, 0, 0), radians(90.0f));
        lbm.voxelize_stl(get_exe_path() + "/stl/NACA-Wide.stl", center, rotation, size, TYPE_S); // set new to Solid

                //reset stuff
        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.x[n] = 0.0f;
                lbm.u.y[n] = 0.0f;
                lbm.u.z[n] = 0.0f;

                lbm.F.x[n] = 0.0f;
                lbm.F.y[n] = 0.0f;
                lbm.F.z[n] = 0.0f;

                lbm.rho[n] = 0.0f;
            }
        }

        lbm.F.reset();
        lbm.rho.write_to_device();
        lbm.u.write_to_device();
        lbm.F.write_to_device();
        lbm.flags.write_to_device();
            lbm.run();// now run it again

And this is what it looks like right after the Revoxelisation, The AoA increased by 5°, the holes at the top right and bottom left are where the airfoil used to be.

And here you can see the white boundary forces that are in the inside of the airfoil:

I get the same issue when moving the airfoil, all the nodes that were not solid at first, have ridiculous ammounts of drag. Has anyone tried revoxelising and getting the forces? And what am I not resetting correctly? I hope everyone was able to enjoy Christmas, Cheers Marius

[Maere05]

I also tested the TIE-fighter from the examples in setup.cpp I get less weird results whilst its turning, since it's only a hull not a solid. but when you look at the drawn lines its still quite weird to me. Before Rotation:

After:

[trparry]

An instantaneous change of 5 deg is quite large. This means that many fluid nodes instantaneously become solid and many solid nodes instantaneously become fluid. So, it isn't surprising that there are large spikes in the forces due to the mid grid bounce back boundary conditions.

If you're going to move a solid boundary you should rotate maybe 0.01 deg per time step or something like that, so you might have to use a small timestep depending on how fast you want to rotate.

Also, I think there is still a challenge to moving boundaries that I haven't had enough time to work on. This is mentioned in #39 and #40 where we have to set the instantaneous velocity of each type_s node according to v = r*omega where omega is the angular velocity and r is the distance to the axis of rotation. Right now you are setting the instantaneous velocity to 0 which isn't physically correct, the fluid that comes in contact with the solid will take on the solid's instantaneous velocity at the point of contact.

I haven't had enough free time to figure this out but I have a feeling its not terribly difficult, you need to iterate over the solid and calculate r for each node, and omega is given by the user.

[Maere05]

I undestand that 5° is ridiculous, it even leaves behind a vaacum. But that was just so it's clearly visible, that once it pitches up there is a ridiculous ammount of force generated from inside of the airfoil. So I think there is something weird with the force calculations, the majority of the forces come from the inside of the solid model, there should not be any force generated from a solid voxel surrounded by solid voxels. something similar happens with intersecting geometries: The mesh of these 2 cubes intersect, so there is no fluid inside of it, therefore the cube shouldnt generate any force from inside of the intersection.

[trparry]

@Maere05 you should look at the latest release of FluidX3D. Improvements were made on the revoxelization which will perhaps make the force vector field make more sense when objects overlap. Also, velocity initialization was added such that the fluid nodes adjacent to a solid take on the instantaneous velocity of the solid during rotation or translation. So, maybe try the overlapping cube test again with the latest release!

Also, it should be noted that the older releases you were using did use a revoxelization algorithm that only voxelized the solid shell, such that there was fluid inside the shell which could explain what you are seeing.

[Maere05]

I already repeated all the tests, the new algorithm is amazingly fast! but it seems to me that there is still something wrong or that I'm missing about this. I don't think there was any fluid inside the Cube in the last test I did, lbm.voxelize_stl(...); was completely solid and voxelize_stl_hull(...); was for the shell. I did most tests with both, the image above was done with lbm.voxelize_stl(...);.

The intersecting Cubes seem to still have forces inside when using lbm.voxelize_mesh_on_device(mesh, TYPE_S); (on Release v2.2) And if I rotate the cube by a degree I get this...

Have you or anyone else changed something in a simulation, like a rotation or translation, and then read out forces? I just cant get that to work somehow and don't really know what im doing wrong.

[trparry]

I have rotated and translated bodies and read forces. I did add a "dwelling period" where I allowed the simulation to run for 10000 steps before moving anything. Also the amount of rotation and translation per timestep was very small (like 0.001 deg per step). So, I don't think I saw spikes in forces like you are seeing. However, the fact that forces are being generated within overlapping regions is a bit strange.

I actually connected matlab simulink to a 2 dof airfoil, where point O1 is a fixed rotational joint and O2 is free to rotate around O1. There is an actuator at O2 that changes the angle between the two bodies. I placed a sine wave on this actuator to see what happens. The results seem to make sense.

https://user-images.githubusercontent.com/41972759/218999086-8efc508d-c68e-4f8e-ae93-047d5671cff0.mp4

[Maere05]

Could you share the code from your main_setup() function? Whenever I rotate something it completely breaks down, even with 10k steps in advance and rotating 0.01° per step. I must be doing something wrong, it quite clearly can work. Thanks in advance

[trparry]

@Maere05 I looked at the way you are calculating forces and I do the same thing:

            //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
    F_d_sim = units.si_F(lbm_force.x); // force in Newton

    const float C_d_sim = F_d_sim / (0.5 *si_rho * A_x * (si_u * si_u));

    //print force components

    print_info("drag force = " + to_string(F_d_sim) + " N");
    fout << to_string(F_d_sim) << std::endl; //send to txt file`

I think this is more a problem of interpretation. People are used to seeing time averaged results on force plots. In OpenFoam we can compare a steady state simulation with a transient one by time averaging the transient results. For example, in a transient sim in OpenFoam we add up all the values over a finite interval such as 50 steps (or whatever interval you've chosen to report results), and then divide this by the total time elapsed in these 50 steps to get the time averaged solution. This removes much of the inherent "noise" that is associated with Karman Vortex streets off of a cylinder for example. While some of this variation is probably legitimate error, a large portion of it is due to fluctuations in forces from the vortex shedding. We just remove this noise so that we can see nice plots, observe the force converging to a steady value, etc. I plotted a time averaged force in FluidX3D and it did indeed look a lot smoother (its basically a low pass filter). This obviously all assumes you are in the turbulent regime.

So, I think whenever you move something, do a very small movement and wait for the time averaged solution to stabilize.

[ProjectPhysX]

The root cause of the large forces after re-voxelization is the lack of initialization of fluid cells that have been solid before. These cells might get a velocity manually assigned for initialization, but the LBM distribution functions in these cells are at an undefined state initially after conversion. Forces are computed directly from the distribution functions.

In version v2.2 with lbm.unvoxelize_mesh_on_device(mesh) I have fixed this, by equilibrium-initializing the distribution functions after solid->fluid conversion. Still re-voxelization steps must be small, less than 1 cell of surface movement per step.

[Maere05]

I did some more testing, with plenty of time between reading the forces and with lbm.unvoxelize_mesh_on_device(mesh) the force Value is correct, the visualization still shows ridiculously long arrows from places that can't have a force on it. Not really an issue for me, since the force vector itself is fine. Some STLs don't voxelize properly for me, mostly in Self-intersecting areas (easily solved with a good 3d modeling program). With good geometry ive had a few cases where a tiny "channel" forms, just a hole of fluid voxels inside the geometry. but these went away with a high enough Resolution.

Thank you for all your work!