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

Moving boundaries clarification #39

Closed trparry closed 1 year ago

trparry commented 1 year ago

Hello, I assume the only difference between the stationary (Type S) boundaries and the moving boundaries is the fact that for the moving boundary case the velocity of the fluid at the boundary nodes are unequal to zero.

So, if I move a "stationary" boundary by revoxelizing the mesh this is actually representing a "slip" condition, because the absolute velocity of the fluid is zero, but the solid has a nonzero absolute velocity. So, the relative velocity between fluid and solid is nonzero (slip).

On the other hand, if I move a "moving boundary" by revoxelizing the mesh, and set the velocity of the fluid boundary nodes equal to the velocity of the solid boundary, then we would represent the "no slip" condition. This is because the relative velocity between fluid and solid is zero, even though the absolute velocity of both is nonzero.

Is this all correct?

Also in setup.cpp you don't have an example of a moving boundary implementation I don't think, you only move stationary boundaries. You recently posted a video of a moving boundary setup on the F1 car by turning the wheels, can you share a code snippet of how you write the solid velocity in the fluid boundary nodes (after revoxelization)? I'm not sure how I would write the solid velocity to the fluid, especially for a rotating object that doesn't have axial symmetry (like a wheel). After this, I think you would call lbm.update_moving_boundaries.

ProjectPhysX commented 1 year ago

Hi @trparry,

Just setting the fluid velocity next to stationary no-slip boundaries to non-zero will not work, as in the next time step the no-slip condition at the stationary boundary almost zeroes out the velocity again.

For moving objects like helicopter blades, you have to enable MOVING_BOUNDARIES and repeatingly:

  1. re-voxelize
  2. set velocity in re-voxelized boundary nodes (lbm.u.read_from_device();, modify boundary node velocity, lbm.u.write_to_device();)
  3. call lbm.update_moving_boundaries() to re-initialize TYPE_MS flags of fluid nodes next to boundary nodes with non-zero velocity
  4. compute a few LBM time steps

Here is the updated F1 setup with ratating wheels. In this case I don't need to re-voxelize, as the wheels stay in place and are almost axially symmetrical and their surface velocity is tangential everywhere.

void main_setup() { // F1 car
    // ######################################################### define simulation box size, viscosity and volume force ############################################################################
    const uint L = 256u; // 2152u on 8x MI200
    const float kmh = 100.0f;
    const float si_u = kmh/3.6f;
    const float si_x = 2.0f;
    const float si_rho = 1.225f;
    const float si_nu = 1.48E-5f;
    const float Re = units.si_Re(si_x, si_u, si_nu);
    print_info("Re = "+to_string(Re));
    const float u = 0.08f;
    const float size = 1.6f*(float)L;
    units.set_m_kg_s(size*2.0f/5.5f, u, 1.0f, si_x, si_u, si_rho);
    const float nu = units.nu(si_nu);
    print_info("1s = "+to_string(units.t(1.0f)));
    LBM lbm(L, L*2u, L/2u, nu);
    // #############################################################################################################################################################################################
    const float3 center = float3(lbm.center().x, 0.525f*size, 0.116f*size);
    lbm.voxelize_stl(get_exe_path()+"../stl/Ferrari_SF71H_V5.stl", center, size); // https://www.thingiverse.com/thing:2990512/files
    const ulong N=lbm.get_N(); uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
        // ########################################################################### define geometry #############################################################################################
        //if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = u;
        if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n] = TYPE_E;
        if(z==0u) lbm.flags[n] = TYPE_S;

        const float3 p = lbm.position(x, y, z);
        const float W = 1.05f*(0.312465f-0.179692f)*(float)Nx;
        const float R = 1.05f*0.5f*(0.361372f-0.255851f)*(float)Ny;
        const float3 FL = float3(0.247597f*(float)Nx, -0.308710f*(float)Ny, -0.260423f*(float)Nz);
        const float3 HL = float3(0.224739f*(float)Nx, 0.210758f*(float)Ny, -0.264461f*(float)Nz);
        const float3 FR = float3(-FL.x, FL.y, FL.z);
        const float3 HR = float3(-HL.x, HL.y, HL.z);
        if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+FL, float3(W, 0.0f, 0.0f), R)) {
            const float3 uW = u/R*float3(0.0f, FL.z-p.z, p.y-FL.y);
            lbm.u.y[n] = uW.y;
            lbm.u.z[n] = uW.z;
        }
        if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+HL, float3(W, 0.0f, 0.0f), R)) {
            const float3 uW = u/R*float3(0.0f, HL.z-p.z, p.y-HL.y);
            lbm.u.y[n] = uW.y;
            lbm.u.z[n] = uW.z;
        }
        if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+FR, float3(W, 0.0f, 0.0f), R)) {
            const float3 uW = u/R*float3(0.0f, FR.z-p.z, p.y-FR.y);
            lbm.u.y[n] = uW.y;
            lbm.u.z[n] = uW.z;
        }
        if((lbm.flags[n]&TYPE_S) && cylinder(x, y, z, lbm.center()+HR, float3(W, 0.0f, 0.0f), R)) {
            const float3 uW = u/R*float3(0.0f, HR.z-p.z, p.y-HR.y);
            lbm.u.y[n] = uW.y;
            lbm.u.z[n] = uW.z;
        }

    }   // #########################################################################################################################################################################################
    key_4 = true;
    //Clock clock;
    //lbm.run(0u);
    //while(lbm.get_t()<=units.t(1.0f)) {
    //  lbm.graphics.set_camera_free(float3(0.779346f*(float)Nx, -0.315650f*(float)Ny, 0.329444f*(float)Nz), -27.0f, 19.0f, 100.0f);
    //  lbm.graphics.write_frame_png(get_exe_path()+"export/a/");
    //  lbm.graphics.set_camera_free(float3(0.556877f*(float)Nx, 0.228191f*(float)Ny, 1.159613f*(float)Nz), 19.0f, 53.0f, 100.0f);
    //  lbm.graphics.write_frame_png(get_exe_path()+"export/b/");
    //  lbm.graphics.set_camera_free(float3(0.220650f*(float)Nx, -0.589529f*(float)Ny, 0.085407f*(float)Nz), -72.0f, 21.0f, 86.0f);
    //  lbm.graphics.write_frame_png(get_exe_path()+"export/c/");
    //  lbm.run(units.t(0.5f/600.0f)); // run LBM in parallel while CPU is voxelizing the next frame
    //}
    //write_file(get_exe_path()+"time.txt", print_time(clock.stop()));
    lbm.run();
} /**/
trparry commented 1 year ago

@ProjectPhysX, ok thanks for the insight! Have you already done a moving boundary with a geometry that requires revoxelization? It would be very helpful to see!

Correct me if I'm wrong: If we were to use a moving wall boundary condition on the tie fighter example (a more complex rotating geometry) this requires revoxelization and a calculation for the velocity of each node of the solid surface. This velocity is then written to the fluid as an instantaneous velocity. In this case every surface node has the same angular velocity around the axis of rotation. The linear velocity of each surface node is v = r * omega where r is the distance from the axis of rotation and omega is the angular velocity. Omega is easy to calculate, but r is a bit more complex based on looking at the setups. Any hints on how to calculate this?

A moving boundary seems to be a more accurate representation to reality than revoxelizing a stationary boundary because we assume that the fluid that comes in contact with the solid will take on the solid's instantaneous velocity at the point of contact. This would also effect the forces on the solid boundary.

ProjectPhysX commented 1 year ago

@trparry I haven't yet tested revoxelization with moving boundaries. give me some time. You're right with your proposal on the tie fighter. It also works to set the velocity to all boundary nodes (also the ones on the inside, might be easier) rather than only at boundary nodes adjacent to fluid nodes. You can calculate float3 r = lbm.position(x, y, z)+lbm.center()-mesh.get_center(); for any solid node at (x, y, z).