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.77k stars 300 forks source link

Issue with Force Vector Positions in FluidX3D Simulation #138

Closed MingLocke closed 4 months ago

MingLocke commented 8 months ago

Dear Dr. Moritz Lehmann,

I hope this message finds you well. I am currently working on simulating the deformation of an underwater wing using FluidX3D. The deformation shapes for the wing within the time range of 0-0.3 seconds are stored in a series of wing_N.stl files (N=0,1,2,...). My objective is to read the new STL file at each time step, calculate and export the position of the force vector on the wing surface, along with the magnitudes of the Fx, Fy, and Fz components, to a CSV file. Below is my setup file:

void main_setup() { // hydrodynamics of a morphing wing;
    // ################################################################## define simulation box size, viscosity and volume force ###################################################################
    const uint3 lbm_N = resolution(float3(2.0f, 1.0f, 1.0f), 2000u); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
    const float si_u = 1.0f;
    const float si_length = 0.4f;
    const float si_T = 0.3f;
    const float si_nu = 1.0E-6f, si_rho = 1000.0f;
    const float lbm_length = 0.8f * (float)lbm_N.y;
    const float lbm_u = 0.1f;
    int current_wing_index = 0;
    units.set_m_kg_s(lbm_length, lbm_u, 1.0f, si_length, si_u, si_rho);
    print_info("Re = " + to_string(to_uint(units.si_Re(si_length, si_u, si_nu))));
    LBM lbm(lbm_N, units.nu(si_nu));
    // ###################################################################################### define geometry ######################################################################################
    Mesh* mesh = read_stl(get_exe_path() + "../stl/stl_series/wing_0.stl");
    mesh->scale(lbm_length / si_length);
    mesh->translate(lbm.center());
    mesh->translate(float3(1.0f - mesh->pmin.x + 0.1f * lbm_length, 0.0f, 1.0f - mesh->pmin.z)); // move mesh forward a bit and to simulation box bottom, keep in mind 1 cell thick box boundaries
    lbm.voxelize_mesh_on_device(mesh, TYPE_S);
    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 (lbm.flags[n] != TYPE_S) lbm.u.x[n] = lbm_u; // initialize x-velocity everywhere except in solid cells
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == Nz - 1u || z == 0u) lbm.flags[n] = TYPE_E; // all other simulation box boundaries are inflow/outflow
        }); 
    lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_Q_CRITERION;
    std::string CSVfilename;
    lbm.graphics.set_camera_centered(-40.0f, 20.0f, 78.0f, 1.25f);

    // ####################################################################### run simulation, export images and data ##########################################################################
    lbm.run(0u); // initialize simulation

    while (lbm.get_t() <= units.t(si_T)) { // main simulation loop

        // Increment the wing index for the next iteration
        current_wing_index++;

        // unvoxelize old stl mesh
        lbm.unvoxelize_mesh_on_device(mesh);

        // load new stl mesh
        std::string wing_file_path = get_exe_path() + "../stl/stl_series/wing_" + std::to_string(current_wing_index) + ".stl";
        Mesh* mesh = read_stl(wing_file_path);
        mesh->scale(lbm_length / si_length);
        mesh->translate(lbm.center());
        mesh->translate(float3(1.0f - mesh->pmin.x + 0.1f * lbm_length, 0.0f, 1.0f - mesh->pmin.z)); // move mesh forward a bit and to simulation box bottom, keep in mind 1 cell thick box boundaries
        lbm.voxelize_mesh_on_device(mesh, TYPE_S);

        // run simulation for 0.001 secs
        lbm.run(units.t(0.001f));
        lbm.graphics.write_frame();

        // export force vector's position and value
        std::cout << "\nCurrent time step: " << units.si_t(lbm.get_t()) << std::endl;
        CSVfilename = get_exe_path() + "lbm_res/node_" + std::to_string(units.si_t(lbm.get_t())) + ".csv";
        std::cout << "Constructed filename: " << CSVfilename << std::endl;

        std::ofstream myoutfile(CSVfilename);

        if (!myoutfile.is_open()) {
            std::cerr << "Unable to open file: " << CSVfilename << std::endl;
            break;
        }

        myoutfile << "X, Y, Z, Fx, Fy, Fz" << std::endl;

        lbm.calculate_force_on_boundaries();
        lbm.F.read_from_device();

        for (ulong n = 0; n < lbm.get_N(); ++n) {
            if (lbm.flags[n] == TYPE_S) {
                const float3 position = lbm.position(n);
                const float3 force = { units.si_F(lbm.F.x[n]), units.si_F(lbm.F.y[n]), units.si_F(lbm.F.z[n]) };

                myoutfile << position.x << ", " << position.y << ", " << position.z << ", ";
                myoutfile << force.x << ", " << force.y << ", " << force.z << std::endl;
            }
        }

        myoutfile.close();
    }
}

The issue I am facing is that the force vectors exported to the CSV file seem to remain at the initial positions and do not update to the deformed locations. To provide a better understanding, I have attached two files:

Attachment 1: A GIF file illustrating the simulation images. Attachment 2: A vector plot generated from the exported CSV file, showcasing the force vectors.

I would greatly appreciate it if you could review my configuration settings in light of the provided attachments and advise on any potential issues. Thank you very much for your valuable assistance.

Best regards

output ForceCSV

ProjectPhysX commented 4 months ago

Hi @MingLocke,

the problem here is that the simulation gets unstable and blows up: in the first GIF about 4 seconds in, it gets fuzzy on the wing edge and then quickly propagates downward. Afterwards, the entire simulation box will be filled with NaNs and you won't get any usable data out. You're at Reynolds number 400k. Please check that you've enabled the SUBGRID extension; but even with that it may get unstable at this large Re and low resolution.

Kind regards, Moritz